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1.0 INTRODUCTION 


NASA's advanced engine programs are aimed at progressively higher efficiencies, 
greater reliability, and longer life. Recent studies have indicated that significant 
engine performance advantages can be achieved by employing advanced seals, and 
dramatic life extensions can also be achieved. Advanced seals are not only required 
to control leakage, but are necessary to control lubricant and coolant flow, prevent 
entrance of contamination, inhibit the mixture of incompatible fluids, and assist in 
the control of rotor response. 

Recognizing the importance and need of advanced seals, NASA, in 1990, embarked 
on a five-year program (Contract NAS3-25644) to provide the U.S. aerospace 
industry with computer codes that would facilitate configuration selection and the 
design and application of advanced seals. 

The program included four principal activities: 

1. Development of a scientific code called SCISEAL, which is a 
Computational Fluid Dynamics (CFD) code capable of producing full 
three-dimensional flow field information for a variety of cylindrical 
configurations. The code is used to enhance understanding of flow 
phenomena and mechanisms, to predict performance of complex 
situations, and to furnish accuracy standards for the industrial codes. 
The SCISEAL code also has the unique capability to produce stiffness 
and damping coefficients that are necessary for rotordynamic 
computations. 

2. Generation of industrial codes for expeditious analysis, design, and 
optimization of turbomachinery seals. The industrial codes consist of a 
series of separate stand-alone codes that were integrated by a 
Knowledge-Based System (KBS). 

3. Production of a KBS that couples the industrial codes with a user 
friendly Graphical User Interface (GUI) that can in the future be 
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integrated with an expert system to assist in seal selection and data 
interpretation and provide design guidance. 

4. Technology transfer via four multiday workshops at NASA facilities 
where the results of the program were presented and information 
exchanged among suppliers and users of advanced seals. A Peer Panel 
also met at the workshops to provide guidance and suggestions to the 
program. 

This final report has been divided into separate volumes, as follows: 


Volume 1: 
Volume 2: 
Volume 3: 
Volume 4: 
Volume 5: 

Volume 6: 


Executive Summary and Description of Knowledge-Based System 
Description of Gas Seal Codes GCYLT and GFACE 
Description of Spiral-Groove Codes SPIRALG and SPIRALI 
Description of Incompressible Seal Codes ICYL and IF ACE 
Description of Seal Dynamics Code DYSEAL and Labyrinth Seal 
Code KTK 

Description of Scientific CFD Code SCISEAL 


This report summarizes the work performed on the Scientific Seals code (SCISEAL) 
under NASA Contract NAS3-25634, during the period of February 1990 to May 1995. 
The work for the first three years (February 1990 - September 1993) consisted of 
development of the single domain version of SCISEAL. Starting from REFLEQS-3D, 
a 3D code developed by CFDRC, several capabilities were added, which included 
Colocated grids, rotating frames, moving grids, high-order differencing and 2-layer 
turbulence models. Two rotordynamic coefficient calculation modules were added 
(whirling rotor method and N-S perturbation method), as well as other seal specific 
modules and boundary conditions. The preprocessor SCIPRE was modified to 
include automated grid generation capabilities for cylindrical seals. During the 
period of October 1994 - May 1995, upon recommendations from CFDRC, MTI and 
the Peer Review Committee, the code was enhanced by incorporation of multi- 
domain capability, ability to treat 2D and 3D problems, and extension of the 
rotordynamics modules to multi-domain format; this is the current status of 
SCISEAL. 
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2.0 CAPABILITIES 


The objective of this task is to develop a 3-D CFD code for the analysis of fluid flows 
and forces in a variety of seals including cylindrical seals (tapered, annular, stepped) 
labyrinth seals, rim seals, tip seals and face seals. This code is to serve as a tool for 
detailed, and accurate analysis of flows in existing seal design as well as new concept 
seals. It can also be used as an accuracy check for the simplified Industrial codes 
which use simplified models for fast turnarounds. Finally, SCISEAL can be used to 
provide detailed flow analyses in secondary flow systems where seals are usually 
coupled with other flow elements such as disk cavities. Such systems need to be 
solved together and in a coupled manner, which is currently beyond the capabilities 
of the seals codes based on the simplified flow models. 

The solver module consists of two separate codes: 

a. The preprocessor SCIPRE, which reads in an input file in a text format 
and converts it into another file that is readable by the main flow 
solver. The preprocessor has capabilities of grid generation, problem 
setup, boundary condition setup and checks for internal consistency for 
the problem as defined by the user in the input file. Detailed 
descriptions of problem setup, command structure etc . is given in the 
Users' Manual for the SCISEAL code. 1 

b. Flow solver SCISEAL is the module which reads in the data files 
created by the preprocessor and does the flow computations. At the 
end of execution of SCISEAL, a number of output files are created that 
can be used for restarting a continuation run and plotting with 
graphical packages; these files also contain information on code 
convergence and integrated seal quantities. A description of the 
various files also is given in the SCISEAL Users' Manual. 

2.1 Salient Features of SCISEAL 

The flow solver, SCISEAL, was written in ANSI FORTRAN-77 with emphasis on 
portability and modularity. The basic flow analysis methodology utilizes advanced 
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numerical techniques for accuracy, efficiency and robustness. Features of the code 
include: 

• Finite volume discretization for integration of Favre-averaged Navier- 
Stokes (N-S) flow equations; 

• Implicit multi-domain treatment with one-to-one and one-to-many 
cell connections at the interfaces; 

• Cartesian, polar, and non-orthogonal Body-Fitted Coordinates (BFC). 

• Colocated (non-staggered) grid; 

• Strong conservative form of momentum equations with Cartesian 
components as dependent variables; 

• Stationary as well as rotating frames of reference for rotary flow 
problems; 

• Pressure-based solution algorithm including a variant of SIMPLEC and 
PISO, which allows the treatment of both incompressible and 
compressible flows; 

• Concentrated and distributed porosity-resistivity technique for 
treatment of internal solid objects; 

• High order spatial differencing schemes (including upwind, central, 
MUSCL and Osher-Chakravarthy) and temporal schemes (Euler 
backwards and Crank-Nicholson); 

• Steady-State and time-accurate solution capability; 

• A space-conserving moving grid formulation that allows the treatment 
of moving-deforming grid systems encountered with whirling rotor; 
and 

• Symmetric whole-field equation solvers based on Stone's implicit 
methods and a conjugate gradient squared solver for linear equation 
solutions. 

The SCISEAL code also has a variety of physical models that are needed in the 
solutions of the flows encountered in seals. These models include: 

• JANNAF property tables for selected species, useful in tracer gas 
simulations with passive scaler transport; 

• Variable viscosity with Sutherland's Law; 

• Advanced turbulence models: 
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- Mixing length model (Baldwin and Lomax) 

- Standard k-e model with wall functions (Launder and Spalding) 

- Low-Reynolds number k-e model (Chien) 

- 2-layer k-e model for rotating flows and for narrow flow passages 

• Isotropic surface roughness treatment; and 

• Comprehensive set of boundary conditions including seal specific 
conditions such as: 

- pre-swirl specifications 

- entrance loss factor. 

22 . Seal Rotordvnamics 

The flow changes in seals due to rotor motion can generate significant fluid force 
changes and affect the overall stability of the rotor system. This effect is introduced 
through the seal rotordynamic coefficients, and two methods are available in 
SCISEAL to calculate the coefficients: 

a. Whirling rotor method: This method uses full CFD solutions for a 
nominally centered rotor whirling in a circular orbit, the solution 
method provides the skew-symmetric set of rotordynamic coefficients 
associated with a centered rotor. 

b. Small perturbation method: The N-S equations are perturbed to 

generate 1st order flow equations that describe the flow changes due to 
small rotor motions. These equations are solved to generate the full 
set of rotordynamic coefficients. This method can be used to treat both 
centered and eccentric seals. 

Both of the methods are fully automated and can be invoked with simple 
commands in the input file. 

Several additional features that are available for seals problems are: 

• Easy grid generation setup for cylindrical seals; 

• Calculation of seal loads, torque and power; and 
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• Facility to read externally created seal grids and use into flow as well as 
rotordynamics modules. 

The following section describes the treatment of flow equations, discretization 
methods and boundary conditions treatment. Descriptions of the physical models 
used in the code and the rotordynamics models follow the numerical models. 

An extensive series of validation and demonstration problems were solved using 
SCISEAL to assess the accuracy of the numerical and physical models and a list of 
these problems is included in Section 4, with a brief description and relevant results 
for each of these problems. 


NAS A/CR— 2004-2 1 3 1 99/VOL6 


6 



3.0 THEORY 


3.1 Geometry and Flow-Domain Modeling 

The SCISEAL code uses a structured grid approach to discretize a given flow 
domain. Several powerful concepts have been built into the code to simplify the 
grid generation procedure, optimize the cell numbers and allow for simpler 
problem definition, these capabilities are: 

a. Generalized non-orthogonal Body-fitted-coordinate (BFC) grids; 

b. Multi-domain capability; and 

c. Internal blockage concept. 

Given a complex flow domain, one or more of these capabilities can be utilized to 
simplify the grid generation process for optimum and accurate flow solutions. This 
section deals with geometry definition, BFC grids, multi-domain approach and the 
internal blockage concept. 

3.1.1 Geometry and Grid Generation 

The seals code uses a finite volume approach, where the flow domain is discretized 
into a number of cells or finite volumes and the flow equations are numerically 
integrated over each finite volume. The discrete representation of the flow domain 
in the computational grid. Furthermore, the grid needs to be of a structured form. 
A grid is considered to be structured if there exists three grid lines (for a 3D problem) 
to identify three distinctive directions and any face of a control volume is on two 
grid lines. In other words, a single (i,j,k) index can be used to identify a cell or a 
point of a structured grid. 

The structured grids can further be divided into single block or multi-block types. 
The single block (or domain) grid has only one set of i, j and k index axes over the 
complete domain. A multi-block grid has several blocks with individual definitions 
of i, j, and k axes, with interfaces at places where the blocks meet. The orientation of 
the axes in the individual blocks is independent of the neighboring blocks. This 
concept is illustrated in Figure 1. 
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Figure 1. Single and Multiblock Grid Concept 

The ability to divide a complex flow problem into several domains is very powerful, 
since this capability allows the user to use grids only on flow areas for optimum cell 
numbers. In addition, the division of flow problem can be made such that the grids 
stay as near orthogonal as possible, i.e. the grid quality can be improved. The 
current capabilities include a one-to-one match as well as a one-to-multi match of 
cells at the interface. This concept is shown in Figure 2. 
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(a) one-to-one interface (b) one-to-many interface 

Figure 2. Types of Interfaces Allowed in SCISEAL 

The ability to locally refine the grids at domain interfaces further enhances the code 
capabilities. Several seal applications have flow domains where a relatively large 
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domain interfaces with a narrow domain, and the multi-to-one interface can be 
used to provide sufficient resolution in the narrow passage, and also reduce the 
number of cells in the large passage through a multi-to-one connection at the 
interface. 

3.1.2 Coordinate Systems 

When solving complex flow problems, definition of curved boundary surfaces with 
Cartesian rectangular grids is difficult and will result in accuracy loss. To be able to 
follow curved lines /surfaces in the flow domain, the so-called body-fitted coordinate 
(BFC) system is necessary. The BFC grids also have 2/3 coordinate axes, however 
they may not necessarily be parallel to any Cartesian axes and may not be locally 
orthogonal to each other. The SCISEAL code has the options to treat Cartesian, 
orthogonal BFC as well as non-orthogonal BFC grids, and such grids can be used in 
conjunction with the multi-domain capability. Cartesian grids require least storage 
and yield highest accuracy, and non-orthogonal BFC grids require the highest 
storage. 

Mathematically, a BFC system can be viewed as a coordinate transformation from 
physical domain to computational domain as illustrated in Figure 3 for a seal sector 
grid. The q and £ coordinates run along the radial and circumferential direction and 
have a direct correlation with the j and k indices. The axial grid lines run along the 
2; direction and are related to the i index. The grid is always represented in a 
Cartesian system in the physical domain. The coordinate transformation converts 
this curvilinear grid to an orthogonal grid in the transformed 

I = 4(x, y, 2 ) , -n = -n [x, y, z) , C = f (x, y, z) 

or x = X{4,T\,Q,y=Y{4,T\,Q,z = Z(4,n, 


space. 


Q 


(i) 
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Figure 3. 



(a) physical domain 



(b) computational domain 

An Illustration of the Transformation from Physical to Computational 
Domain for a Seal Sector Grid 
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Therefore, each point on a physical domain can be identified by a triplet (x, y, z). 
This triplet is a function of t\ and £ or is associated with a particular (i, j, k) index. 
Thus SCISEAL works on a BFC grid system expressed on a base Cartesian coordinate 
system . Note that a physical domain is always transformed to a rectangular domain 
as shown in Figure 3. 

For future reference, a brief introduction of the BFC coordinate system is given 
below. A comprehensive description of BFC coordinate system can be found in the 
work of Thompson et al,2 For derivation convenience, (£1 ^2/ ^ 3 ) is used to replace 
(5, r\, Q. In a BFC coordinate system, (§ 1 , %i), the covariant base vectors are 

defined as 


Zj = isT.'i - 1 ' 2 ' 3 




(2) 


♦ 44 

where ^ is the displacement vector and is equal to x i + y j + z k , and contra variant 
base vectors are defined as 


5*1 _ g 2 x g 3 . *2 _ S 3 x . *3 _ jjl x g 2 

C— J ,C— J , C — J 


( 3 ) 


where J is the Jacobian defined as J = (2 3 x § 2 ) * £3 • Basically, £y is a base vector along 
coordinate line and V is a base vector normal to the surface formed by J. and 
£ k (i*j * k ) . But note that £y and are not unit bases and unit bases can be defined 
as 



2 ; 




( 4 ) 


and 


NAS A/CR— 2004-2 1 3 1 99/VOL6 


11 



(no summation on j) 


(5) 



where hj is usually called scale factors. 

It is easy to show that the covariant and contravariant base vectors satisfy the basic 
relationship: 


• V = Sij ( 6 ) 

where is the Kronecker delta. 

In the BFC system, the gradient, divergence, curl, and Laplacian operators can be 
expressed in conservative form as 


Gradient: 




(7) 


Divergence : 


Curl: 




VxV = 4- d 




(le k v) 

- k (ji k xv) 


(8) 

(9) 


Laplacia n: V 2 f = I f ? ■ tf\ (10) 

The geometric meaning of the above quantities can be easily explained on 2D 
curvilinear coordinates at point P as shown in Figure 4. For a BFC grid, the Jacobian 
is simply the volume of each cell. 
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Figure 4. Geometrical Meaning of Covariant and Contravariant Bases for a 2D 
BFC System 

3.1.3 Blockage Concept 

In many engineering problems the boundaries of the flow problems can be very 
complicated, and there can be internal solid obstacles. In some cases such obstacles 
also need to be included for performing conjugate heat transfer analysis. Treatment 
of such internal objects can either be done through using the multi-domain 
approach or through using the so-called blockage concept. With blockage capability, 
several parts of the flow domain can be excluded from the calculations or such 
parts can be treated as solids, to be included in the flow calculations for conjugate 
heat transfer. The blocked region concept is particularly useful in the treatment of 
seals such as stepped cylindrical seals and labyrinth seals. A labyrinth seal is shown 
in a section along the axis in Figure 5. 
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4115/412-3 


Figure 5. Cross-Sectional View of a Labyrinth Seal with Blocked Cells for Teeth 
Region 

In such a problem, a BFC grid can also be employed to cover the computational 
domain, however, such a grid is difficult to generate, and will have very skewed 
cells, especially near the comers, which can affect the solution accuracy. With the 
blocked-cell concept a Cartesian orthogonal grid can be utilized. The computational 
cells in the labyrinth seal teeth now can be treated as blocked-off or included in the 
energy equation for conjugate heat transfer. In this case, use of orthogonal grid both 
yields ease of grid generation, better accuracy and computational speed, although the 
blocked cells do represent a penalty in terms of computational effort needed to treat 
non-flow cells. 

3.2 Basic Governing Equations 

The fluid flows are simulated by numerically solving partial differential equations 
that govern the transport of flow variables. These variables include mass, 
momentum, energy, turbulence quantities, and mixture fractions. The variables for 
which transport equations have to be solved will depend on the nature of the flow 
problem. In this section, the basic governing equations for the conservation of 
mass, momentum, energy are presented. The user will be introduced to the 
transport equations for other flow variables such as turbulence quantities in a 
following sub-section. 

SCISEAL employs conservative finite-volume methodology and accordingly all the 
governing equations are expressed in conservative form. Cartesian coordinate 
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system and tensor notation are generally employed in which repeated indices imply 
summation over all coordinate directions. The user should note that SCISEAL 
solves the governing equations in cylindrical coordinate system for 2D axisymmetric 
flow problems. Presented here are the flow equations that correspond to the 3D 
representation. 

3.2.1 Continuity Equation 

In any fluid flow in which the continuum hypothesis holds mass conservation can 
be expressed as 


dp d . 

Bt + ^( pu i} = 0 UD 

where uj is the j th Cartesian component of the instantaneous velocity, and p is the 
fluid density. 

3.2.2 Momentum Equations 

These equations are derived from the law of conservation of momentum. 

¥ {pMi ) + ^( pM '“/) = '^ + 'g iTj+pfi (12) 

In the preceding equation p is the static pressure, X\j is the viscous stress tensor and 
fi is the body force. The viscous shear stress, Tfy, is related to the mean shear rate 
(strain rate) tensor, y if/ by 



(13) 


where is the Kronecker delta and 



( 14 ) 
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and /i is the fluid dynamic viscosity which, in general, can be a function of 
temperature, species composition, position, shear rate, etc. If viscosity does not 
depend on shear rate, the fluid is called Newtonian and r i; - is a linear function of 

velocity gradients 



(15) 


Substitution of Equation (15) in Equation (12) results in the Navier-Stokes Equations 



dUj\ 

dXi) 


2 ZuiA 
- 3 


+ Pfi 


(16) 


3.2.3 Energy Equation 

The equation for the conservation of energy can take several forms and different 
forms are suitable for different classes of problems. In SCISEAL, the user may 
choose either static enthalpy or stagnation enthalpy depending on the application. 
The static enthalpy form of the energy equation can be expressed as 


dr i \ d r i \ tyj dp dp dUj 

) + ± (P“ r h ) = ■ ^ + w + u i s: + 


(17) 


Here, qj is the j-component of the heat flux. Fourier's Law is employed to model the 
heat flux 


r E 

q r' K ate. 


(18) 


where K is the thermal conductivity. By algebraic manipulation, the heat flux can 
be expressed in terms of the enthalpy gradient as shown below. 




k_ ar __ a/i_ 

C P rtej- r dx f 


(19) 
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Here C v is the fluid specific heat and T is known as the diffusion coefficient and is 
related to the fluid dynamic viscosity through the Prandtl number a. 


r= 



(7 


The static enthalpy equation can be rewritten as 


( 20 ) 


dt 


C P h ) + 


a , M a 
&i pu > h ) = 


Jdh\ 

dX: 


)i 


dp dp 
+ dt +U ifo'; 


dUi 


+ T; 


V dXj 


( 21 ) 


Note that the above equation is not strictly conservative by its nature and is 
recommended for incompressible and low Mach number flows. On the other hand, 
the total enthalpy form of the energy equation is fully conservative and is 
recommended for high speed compressible flows. The total enthalpy H is defined as 


H = h + (22) 

and the governing equation for H is obtained by adding the fluid kinetic energy 
equation to the static enthalpy equation 




(23) 


3.2.4 Favre Averaged Equations 

The fundamental equations of fluid dynamics that have been introduced in the 
preceding sections are, in general, applicable to Newtonian fluid flow under steady 
or transient, incompressible or compressible, laminar, transitional or turbulent 
conditions. The nonlinearity of the Navier-Stokes equations, coupled with the 
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complexity of the boundary conditions, makes it impossible to obtain analytical 
solutions for all but a limited number of flows of engineering interest Hence one is 
forced to resort to approximate or numerical methods. As most engineering 
applications only require time-mean quantities, the Navier-Stokes equations are 
usually averaged over time or ensemble of statistically equivalent flows to yield 
averaged equations. In the averaging process, a flow quantity 0 is decomposed into 
mean and fluctuating parts. The following two types of averaging are generally 
used. 


Reynolds (or time) Averaging: 

0 = 0 + 0 where 0 ~ (1/T) f 0 dt 
Favre (or density) Averaging: 


^ fw ^ ___ _ 

0=^+0 where 0 = p0/p 


(24) 


(25) 


Note that overbar denotes Reynolds averaging while tilde denotes Favre averaging. 
The time period of averaging, T, should be large compared to the fluctuation time 
scale so that mean quantities are stationary over a number of samples. The user also 
should bear in mind that the mean quantities can vary in time on a scale much 
larger than T. 


Applying the Favre averaging procedure to Equation (11), we get the Favre-averaged 
continuity equation 


dp d 

w + s;H =0 

Similarly when the Navier-Stokes equations (Equation 12) are averaged, we obtain 
the Favre-averaged Navier-Stokes (FANS) equations given below. (For detailed 
derivation, see Cebeci and Smith3). 





A* 


'dUj dUj 2 

< W j + W i '3 dx m 



(27) 
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The FANS equations contain less information than the full NS equations, but have 

additional unknown terms -pupi* called the Reynolds stresses. These correlations 
between the fluctuating components arise in the averaging process, and need to be 
modeled to achieve closure of the FANS equations. All the turbulence models 
available in SCISEAL employ the generalized Boussinesq eddy viscosity concept in 

which the Reynolds stress - p u^u* is treated as a linear function of the mean strain 
rate 




du- dH: 7 1 2 


dx i 


dx t 3 dx m 


(28) 


Here fit is known as the turbulent eddy viscosity and k is half the trace of the 
Reynolds stress tensor. 


k = 



(29) 


By substituting Equation (28) in Equation (27), we obtain the modeled FANS 
equations 


d f - d dp d 


7 L 


du- dUj 2 dUj 


\ 


(' i+ ' , 'k + s:i 


dx i 3 dx m l i 




2 9 riA 

-js;W 


(30) 


When the averaging procedure is applied to the static enthalpy equation, additional 
terms containing fluctuating enthalpy and velocity components appear which are 
generally modeled with the Boussinesq concept as follows 


-P“i h =r ‘&- j 01 ) 

The term r t is known as the turbulent or eddy diffusivity and is related to the 
turbulent viscosity through the turbulent Prandtl number Of. 
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( 32 ) 


The mean total enthalpy equation is derived by adding the mean kinetic energy 
equation to the mean static enthalpy equation. Without any formal derivation, the 
modeled mean enthalpy equations are 


dph d ,-~ 

“ar + a? 




dU: ( dU: dU: 2 ^ U 


(33) 



Various models differ in the way p t is estimated. The models employed in 
SCISEAL are discussed in a subsequent section. 

3.3 Discretization Methods 

The fluid flows are governed by several physical conservation laws and these laws 
can be written as Partial Differential Equations (PDE's) as presented in Section 3.2 A 
numerical method to solve these PDE's consists of the discretization of the PDE's on 
a given grid, formation of corresponding linearized algebraic equations, and the 
solution of the algebraic equations. This way, a final set of discrete numbers on a 
grid is obtained which represents the numerical solution of the PDE's. In this 
section, the discretization of the governing equations is presented. The finite- 
volume approach is adopted due to its attractive capability of preserving the 
conservation property. In the following, the storage locations of the dependent 
variables are first discussed. The discretization process then follows in more detail. 
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3.3.1 Staggered Versus Colocated Grid Approach 

A detailed description of this method was provided in the First Interim Report 4 and 
only a short description is given below for brevity. The staggered grid approach was 
widely used for incompressible flow simulations in the past. With the staggered 
grid approach, proposed by Harlow and WelchS, the velocity components are stored 
at positions between the pressure nodes as illustrated in Figure 6a. Such an 
approach ensures that the pressure is readily available for momentum equations 
and velocity components are available for the continuity equation without 
interpolation. As a result, a proper pressure coupling is guaranteed and the well- 
known checkerboard instability is prevented. However, the disadvantages of the 
staggered grid approach are well known and the following are just two examples. 

• It is not easy to extend the staggered grid approach to non-orthogonal 
curvilinear grids. Several proposed extensions are extremely complex 
to apply and can cause loss of accuracy. 

• Many state-of-the-art CFD methodologies are difficult to apply with a 
staggered grid approach such as multigrid, local grid refinement, and 
multi-zoning. 

The current state-of-the-art approach is based on the Colocated grid (or non- 
stagger ed grid) proposed by Rhie and Chow 6 and later by Peric. 7 The grid 
arrangement of this approach is illustrated in Figure 6b. This Colocated grid 
approach has many advantages over the staggered grid approach and is used in the 
seals code. This approach evaluates the cell face velocity using momentum 
equations and consequently the cell face velocity is directly linked to a third-order 
pressure derivative term. This linkage ensures that the checkerboard instability is 
eliminated. With a Colocated grid approach and Cartesian components as 
dependent velocity variables, many coding complications are avoided and more 
accurate solutions can be obtained, particularly for viscous flows. 
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Figure 6. Illustration of Variable Storage Locations for: (a) Staggered; and 
(b) Colocated Grid 


3.3.2 A General Convection-Diffusion Equation 

It is noted that all governing equations, except continuity equation, can be expressed 
in a general form as 


dp(J) dpucj) dpzx|> dpw;<j> 


dt 


+ 


dx 


+ 


dy + 3z 



(35) 


where <p can stand for Cartesian velocity components, total enthalpy, turbulence 
kinetic energy, mass concentration, etc. r is the effective diffusivity and is the 
source term. Therefore, the above equation, in convection-diffusion form, will be 
considered for discretization. The continuity equation has a different form, and will 
be discussed in a later section. 

First of all. Equation (35) needs to be transformed to BFC coordinates using the new 
independent variable tfx, y , z), r\(x, y, z), and y, z). Without detailed derivation, 
it is sufficient to write down the transformed equation in (£, r\, Q coordinates as 
follows 
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(36) 


with 


$t<j p « + t 4 (/p( ^ 


M A \ 


r Jg> k 


9<J> 


+ s„ 


gft = & . g* 


(37) 


The discretization involves an integration of Equation (36) over a control volume as 
shown in Figure 7. That is: 


<ttt> 




r/« 


;/t 3<t> 




+ s , 


J dl % dr I dC (38) 


In the next several sections, individual terms of the above integration will be 
discussed. 
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3.3.3 Transient Term 

Consider the discretization of the transient term of Equation (36), i.e., 


T - 


y (/p<t>)d|dridf = 


p4>v - 

_ 


( 39 ) 


In the above, superscript “o" stands for a value at an old time, variables without " o " 
superscript are at the new time, and V stands for the volume. Note that the above 
discretization holds true for both Euler first-order and Crank-Nicholson second- 
order schemes. 


3.3.4 Convection Term and Different Convection Schemes 

By defining a physical contravariant velocity component U k such that V = U k & k , the 
convection term can be rewritten as 


c= /i [/p( ^^ 1= 7i[i puJ$ 


+ u 

} dr\ 


hi pU2 \ 


+ 11 . 
/ 


P U 3 4> 


(40) 


Integration of the above term over a control volume gives 


C = C e - C w + C n - C s + C h - Ci - G e ( j) e - G w 4> w + G n bn ' G s <^ s + - G/bi 


(41) 


with G e defined, for example, as 


C,,(±p U ') i (42) 

where G's represent the mass flux through a face of the control volume. The upper- 
case subscripts, W, E, S, N, L and H are used in this report to denote the neighboring 
cell-centers on the west, east, south, north, low, and high sides of the control 
volume shown in Figure 8. The lower-case subscripts w, e, s,n, Z, and h, are used to 
represent the corresponding cell faces of the control volume. The evaluation of the 
mass fluxes, G, will be described in Section 3.4 while the evaluation of <j> at control 
volume faces is described next. 
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For ease of illustration, consider the 2D control volume as shown in Figure 8 
Because the flow variable (p is available only at the cell-centers, the cell-face values 
of <p need to be interpolated. Various interpolation schemes with varying levels of 
numerical accuracy and stability are in use today. In SCISEAL the user has a choice 
of several popular schemes, each of which is illustrated below in the evaluation of 
<p e/ the value at the east cell face. 



Figure 8. A 2D Stencil for the Discretization of Convection and Diffusion Terms 


3.3.4.1 First Order Upwind Scheme First order upwind approach will evaluate <p e 
using either <pp or fa depending on the flow direction at point e. Mathematically it 
can be expressed as 


I G e $ p if G e >0 

y G e 4>E if Gg < 0 (43) 
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or in more convenient form as 


UP fa + fy 

~~ W 2 



4>e ~ 4>p 
~1 


(44) 


This scheme has first-order accuracy and is one of the most stable schemes. 

3.3.4.2 Central Difference Scheme Pure central difference approach will evaluate (j> e 
by averaging the values at point P and E. That is 


r CN 


rg 2 


(45) 


It is widely known that pure central difference scheme may cause unphysical 
oscillations. For most problems, some damping (or artificial viscosity) is needed for 
stability. In practice, central difference scheme with damping is constructed as 


C e = d e C“ p + (l-d e )Ce N (46) 

with d e , the damping coefficient, representing the fraction of upwind scheme used. 
Equation (46) can be re-written as 


C e = G e 


+<I ) E 


“ I G e 




(47) 


Comparing Equation (47) with Equations (45) and (44), it is clear why the scheme is 
called central differencing plus damping. d e = 0 yields the pure central scheme, 
while d e = 1.0 results in the upwind scheme. In this report, central difference 
scheme often refers to pure central plus the damping term as given above. 

3.3.4.3 Second Order Upwind Scheme In this scheme, the cell face value is 
evaluated using two upstream nodes. The cell convective flux is given by 
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(48) 


I G «(f “ 2 ^ w ) tf Ge >0 
\ Gg (f “ I ^Ee) tf G « < ^ 


3.3.4.4 Smart Scheme with Minmod Limiter For the central differencing scheme, 
the damping coefficient is constant for the entire domain. There are many flow 
situations where damping is needed only in certain limited regions. Therefore, 
smart scheme is designed to adaptively calculate damping coefficient depending on 
the local flow property. The minmod limiter can be used to obtain damping 
coefficient, d e as: 


d g = minimode [7, yj 


(49) 


with 

U ~ - U , 1T • (r*\ 

~1e = ipx and U = st g n ( G e) (50) 

The MINMOD fimction is defined as 

MINMOD (a, p) = sign (a) max [o, min 0 a [ p)] 

This scheme will reduce to upwind scheme if there exists local extreme such as 
across shocks (y e < 0), to central difference scheme when y e > 1 , and second order 
upwind scheme for 0 <y e < 1. 

3.3.4.S Other High-Order Schemes For compressible flows with shocks, several 
high-order schemes with limiters are very accurate for shock capturing. Three such 
schemes are introduced in the following. 

Osher-Chakravarthv Scheme: Osher-Chakravarthy scheme evaluates the damping 
coefficient as 


d e = ?—^-minmod{ p, yj + ^-^-minmod (I, pyj 


(51) 
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with t| = j and p = y— 


Roe's Superbee Scheme: Roe's Superbee scheme is used to obtain the damping 
coefficient as 


d e = Max [o, min (2y, 1), min [y, 2)] (52) 

van Leer's MUSCL Scheme: van Leer's MUSCL scheme is used to obtain the 
damping coefficient as 


A lYel+Ye 

e "(3+ Ye) 


(53) 


3.3.5 Diffusion Terms 

The diffusion term in Equation (36) can be split into two parts: main diffusion 
(j = k); and cross diffusion (j * k). Lefs consider main diffusion first, i.e. 


n fc — I 

u M-~ 


J 


r Js 


kk a< t> 




k= 1,2, or 3 


(54) 


Without loss of generality, k = 2 term will be used for derivation. Integration of D M 
term over a control volume, as shown in Figure 7, leads to 


1JT 


Dm 1 dl;dr\d^ = 


r Is 


ud± 
K 


k 1 


J w 


(55) 


It is easy to show that 


ro.il = df. = d 

J hjsin cq 


(56) 
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where A is the area of the control volume face along the § direction and cci 
represents the angle between vector e 1 and the plane formed by e 2 and e 3 . By 
defining 


we have 



r. a 


i 


h x sin « 1 


(57) 


JJJ Dl J dt dr! d( = D] (<P E - - of (#, -*„)=- (of + + D? ^ ^ (58) 

Therefore, D £, the main diffusion coefficient needs to be evaluated at the faces of 
each control volume. 

The cross diffusion term is written as 


jjft - X - d _ 

c iWk 


JTg > * 


d<|) 

*5 


j*k 


(59) 


Consider D 22 for illustration purposes 

1/1 D ? / wf-K £),-K $). « 


By defining D 22 = J Jg 22 and assuming ^ ^ ($p + fe + fe + fes) / etc -> we have 


_ 1 


JJJ J d%dr\d£ = D 22 (fe + fe E - <j) s - <J>se) ' &w (fe + ^nw " fe “ few) (^1) 


Other cross terms can be similarly discretized. 

3.3.6 Source Term 

If the source term is a function of 0 itself, it is linearized as 


NAS A/CR— 2004-2 1 3 1 99/VOL6 


29 



= S F ip + S u 


(62) 


such that S p is negative. In general, both S p and S u can be functions of <p. They are 
evaluated using the latest available value of (f>, which is generally taken to be the 
value of 0 at the end of previous iteration. The linearized source term is integrated 
over the control volume which results in 


S — Sp (pp + Sjj 

(63) 

where the coefficients Sp and Su are given as follows 


S P =VS P 

(64a) 

s u = vs u 

(64b) 


3.3.7 Finite Difference Equation 

In Sections 3.3.3 to 3.3.5, the general convection-diffusion Equation (36) has been 
discretized term by term over a control volume. If the numerically integrated 
transient, convection, diffusion and source terms are assembled together, it results 
in the following linear equation. 


a p <pp = + + a s < h + a b$N + fl iA + fl H^i (65) 

+ fl SW$SW + a SE$SE + fl NW^Ww + fl NE^VE 
+ a LS$LS + ZlnAn + %sfos + fl HN^N 

+ a wrf>WL + fl WH^WH + fl ElAl + a EH$EH 
+ Sjj 

The coefficients ap, aw, etc. are known as link coefficients, and the Equation (65) is 
known as a finite difference equation (FDE). This equation is the discrete equivalent 
of the continuous flow transport equation we started with. To summarize, the PDE 
represented by Equation (36), when numerically integrated using the finite-volume 
methodology, results in the FDE given by Equation (65). 
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This FDE, in general, is nonlinear because the link coefficients themselves are 
functions of <pp, (pw, etc. When an FDE is formulated for each computational cell, it 
results in a set of coupled nonlinear algebraic equations. No direct matrix inversion 
method is available to solve a set of nonlinear algebraic equations. Therefore an 
iterative procedure is employed in SCISEAL at every time step. A linear FDE is 
formed by evaluating the link coefficients with the values of <f> available at the end 
of previous iteration. 


a k P 4 +1 =^ k nb i b +1 + S k u ( 66 ) 

Here, the compact notations a n b and <p n b are used to represent the link coefficients 
and the values of flow variable corresponding to the neighboring grid points. The 
subscripts k and fc+2 denote the previous and current iteration numbers 
respectively. When the linear set (66) is solved, we have an improved estimate for 
(f>. This improved estimate is used to update the link coefficients ap, a n b and Su and 
the linear set is solved again. The iterative procedure is repeated until a converged 
solution is obtained. 

3.3.8 Pressure Gradient Term 

The pressure and velocity fields are strongly coupled in fluid flows. For this and 
other reasons that will become obvious to the user in the next chapter, the pressure 
gradient term appearing in the momentum equations is treated differently from 
other source terms. 

Let us consider the pressure gradient term of the x-momentum equation. In a BFC 
system, it is represented as 


dp_ 

dx 





where r\ x and Cx are given by 


(67) 



Vx 

£ 



( 68 ) 
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Equation (67) may be represented in a compact form as shown below. 


iP (Jc 
dx~ 04 1 


( 69 ) 


When 5p / dx is integrated over the control volume shown in Figure 7, it results in 
the following 


err 

J i)d$dndC 

= a 2x + a 2 * + a 3x Vgj 


(70) 


where, A**, A 2x and A$ x represent projected areas of control volume faces, e. g., A^ x 
is the area of face of the control volume projected normal to the x-direction. 
Vg>, V^p and Vqp represent pressure differences which, when evaluated at the cell 
center, are given by 

v fp=Pe-p„ 

V n p = p n-Ps (71) 

V ip=pk-P( 


where p e , p s etc. are cell-face pressures that need to be interpolated from neighboring 
nodal pressures. Equation (70) may be rewritten in the following compact form 

^%) =AkXV kP ( 72 ) 

Here, 8p/8x is the discrete equivalent of the continuous pressure gradient dp/dx . 
Similar expressions can be obtained for dp/dy and dpjdz. 
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3.4 Discretization of Mass Conservation and Mass Flux Evalu ation 


As mentioned before, mass conservation equation is a special one which can not be 
written as the general convection-diffusion form (Equation (36)). Moreover, in a 
pressure-based method mass conservation is used to determine pressure field. For 
this purpose, the mass conservation equation needs special attention. The general 
mass conservation equation can be transformed into BFC coordinates as 

<73> 


Therefore, integration over a control volume yields 

pV ~a£ — - + C e -G w + G n -G s + C h -G l = 0 


(74) 


With G = (— p U 1 ) for example. Note that G's are the mass flux through a control 
* U 1 ! Je 

volume face as mentioned in Section 3.3.4. 


The next task is to evaluate mass fluxes G's at control volume faces such that the 
checkerboard instability is eliminated for a Colocated variable arrangement. This is 
achieved by averaging momentum equation to the cell faces and relating the cell 
face velocity directly to the local pressure gradient. A brief description of the 
procedure is outlined below. 


The discretized x-momentum equation can be expressed at a cell center P as (see 
Figure 7) 




= 2afeu ni ,-V 

no 


pIS. 


f pV' 

I At K 


U\ 


+ S“ 


(75) 


where nb refers to all the neighboring cells and the transient term is explicitly 
expressed to ensure that the momentum equation at cell face includes the effect of a 
time step. Dividing the above equations by a“ yields 
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Unb ' d P V3z + cd r u °P + a* 


c = £ 


d u = — 
' a P ait 


The above momentum equation is at P-cell but in reality it applies to every point. 
Therefore, at cell face /, for example, we have 

f> - “iy *r -(5 f : “■*) ■ d V li), “f * (f ), ™ 


Since we do not know how to evaluate (a% b ) j , for example, the following quantities 

at / will be obtained by averaging the same quantity from two neighboring cell center 
points. 


ju (y^nb u ) (Su\ 

d f' 


Therefore, by using Equation (76) we have 


P a$ a$\ f a“ 


= (l+c^p + ^^J -cd^ 


where overbar means average from points P and £. A second-order accurate linear 
interpolation procedure is used in SCISEAL to get the various averaged quantities at 
cell face /. Equation (78) can be written as: 


P + «*“) “/ = p + ^>, + *j; ($)-<*;«? - d; + a*; * 
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and the above equation is further reduced to 


with 



( 82 ) 


d u 


i-cd ; 


(83) 


Now it is clear from the above equations that the cell face velocity is obtained from 
an average of the two neighboring point velocity plus a second-order pressure 
correction and a second-order previous time velocity correction. The pressure term 
serves as the mechanism to avoid the checkerboard problem and the previous time 
velocity serves as a mechanism to obtain a time-independent steady solution. 


The y-component and z-component velocities can be similarly evaluated and the 
contravariant velocity component at the cell face can be thus evaluated as 

U} = u r F u + v r Fi y + w f -F u (84) 

with 

F, =h.e! ■ i ; F, = h.e - J ; F, = h, e 1 ■ k (85) 

lx 1 ly 1 J lz 1 


and the grid metrics, Fi x etc. are evaluated at the cell face. 

3.5 Pressure-Correction Equation 

As outlined earlier, each flow variable is governed by a partial differential equation 
(PDE) which is numerically solved to obtain a discrete solution for that variable. 
The three momentum equations yield the three cartesian components of velocity. 
Even though pressure is an important flow variable, no governing PDE for pressure 
is presented. Pressure-based methods utilize the continuity equation to formulate 
an equation for pressure. Two methods to achieve the velocity-pressure coupling 
are available in SCISEAL: 
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• SIMPLEC 

• PISO 

The basic framework of each method is briefly described below. 

3.5.1 SIMPLEC Algorithm 

SIMPLEC stands for "Semi-Implicit Method for Pressure-Linked Equations 
Consistent", and is an enhancement to the well known SIMPLE algorithm^. In 
SIMPLEC, which was originally proposed by Van Doormal and Raithby?, an 
equation for pressure-correction is derived from the continuity equation. 

The finite-difference form of the ^-momentum equation can be written as 

a p u p = {Za nb u nb + V k p p (86) 

where the subscript P refers to cell center P. The pressure field should be provided to 
solve Equation (86) for u . However, the pressure field is not known a priori. If the 
preceding equation is solved with a guessed pressure (or the latest available pressure 
in an iterative procedure) p*, it will yield velocity u* which satisfies the following 
equation. 

Up Up = {£a nb u nb + S u } p - Af v k p‘ p (87) 

In general, u* will not satisfy continuity. The strategy is to find corrections to u* and 
p* so that an improved solution can be obtained. Let u' and p' stand for the 
corrections. 


* * 
u = u +u 

(88) 

* ' 

p= p +p 

(89) 


An expression for 




can be obtained by subtracting Equation (87) from (86). 


a P u p = {Za nb u' ni } p -Af V k p' p 


(90) 
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Because our aim is to find an incrementally better solution, we will approximate u ^ 

/ ' 
by Up giving us an expression for u p 


Up - d k p V k p p 


where 


A _ 

dp~ 


i kx 


a P'^ a nb 


(91) 


(92) 


We can also obtain an expression for u' e , the velocity correction at the east cell face, 
by averaging from the nodal corrections 


' (—£ ' 

u =- 4+4 v.p 

e V P E ) k^e 


(93) 


and the gradient of the pressure correction is evaluated at the cell face e. 
The discretized continuity equation is 

Pp ■ Pp 


At 


+ G e- G w + G n~ G s + G h~ G t - 0 


(94) 


The convective fluxes G e etc. can be represented as 


G = G*+G' (95) 

Where G* represent the fluxes calculated using u*, P* etc., and G' is the correction to 
the fluxes, and is to be evaluated. 

Substituting in Equation (94) and rearranging, the discretized continuity equation is 


Pp V P 
At 


+ G - G„, + G„ - G e + Gl - Gt — S, 


(96) 


Where S m represents the mass correction or mass source in the control volume 
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( 97 ) 


S m = - 

ttl 


0,0 * 


Pp ^p ” Pp ^p * * * * * _* 

— + G e “ G w + G n “ G s + G fe ' G t 


At 


By using Equation (42) in Equation (95) we get an expression for 






(98) 


In the derivation of the preceding equations, products of primes have been 
neglected. 

For incompressible flows p' is zero. For compressible flows p' is estimated as 



An equation of state is used to estimate dp/dp. If the fluid is a perfect gas 

dp i 

dp RT 


(99) 


( 100 ) 


The density correction at the east cell face, p', is estimated from Pp and Pp using the 

i jl 

same scheme used to estimate p e . u e is obtained from the cartesian velocity 
corrections. 


Vi = u e F Tx + v e F ly + W e F lz 


( 101 ) 


where V\ x > & c ' are hie grid metrics at cell face e. Equation (93) is used to obtain u e 

and v e and w e are obtained similarly starting with y- and z-momentum equations 
following the procedure outlined in this section. 

Expressions for the contravariant velocity corrections and fluid density corrections 
at all the cell faces are obtained and substituted in Equation (96) yielding an FDE for 
the pressure correction 


NAS A/CR— 2004-2 1 3 1 99/VOL6 


38 



a p Pp - fl wPw + a E Pe + a S Vs +ci nVn + a L Pl + U hPh+ 


( 102 ) 


The SIMPLEC procedure can be summarized as follows: 

1. Guess a pressure field p*. 

2. Obtain u*, v* , and w* by solving discretized momentum equations such as 
Equation (87). 

3. Evaluate G* from p*, u*, etc. using the procedure outlined in Section 3.4. 

4. Evaluate mass source from Equation (97). 

5. Obtain p by solving Equation (102). 

6. Use p' to correct the pressure and velocity fields using Equations (88) and (89). 

7. Solve the discretized equations for other flow variables including turbulence 
quantities. 

8. Go to step 2 and repeat the procedure until convergence is obtained. 

A variation of the procedure outlined above is employed in SCISEAL. After Step 6, 
instead of executing Step 7, Steps 3, 4 , 5, and 6 can be repeated a few times. That is, 
only the pressure-correction equation is solved a few more times and the pressure, 
velocity and density fields are updated. These intermediate iterations are called the 
'continuity or pressure-correction iterations'. They are found to enhance overall 
convergence for flow of compressible fluids. 

3.5.2 PISO Algorithm 

PISO stands for Pressure-Implicit with Splitting of Operators. PISO employs a non- 
iterative time marching method to handle the velocity-pressure coupling. It is 
essentially a predictor-corrector method in which the velocity, pressure and density 
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fields from the previous time level are incrementally corrected to obtain the flow 
field at a new time level. PISO was developed by Issa^o. 

In order to explain the framework of PISO, it is convenient to recast the discretized 
continuity equation in the following form 

a o 

o V - o V 

— + VAi ( pu; ) = 0 (103) 


where A,- is the finite difference equivalent of the divergence operator and subscript 
"i" indicates cartesian x, y, and z directions. Equation (103) can be rewritten as 



p-p pi V°\ 
At ' At[ V] 


(104) 


The discretized momentum equations are represented as 

pu { V-p° ul 

— h UpUi — £ a n b u i f nb~ VA i P + Sjj (105) 


With some rearrangement of terms, the preceding equation can be written as 

I V a p\ P° u i v ° 

\-^ + jJP u i= I:a nb u i,nb- VA iP + S U + —^— (106) 

Note that both the continuity and momentum equations are numerically integrated 
over the control volume with P as its center. The subscript P on p, p, Ui and V is 
omitted for convenience. 


Let us represent the dependence of p on p as 

p = pf(p,T) (107) 

For an incompressible fluid, / = 0. For a perfect gas. 


f(P'V = 


1 

RT 


(108) 
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PISO algorithm with one predictor and two corrector steps is illustrated below. In 
principle, there is no limit on the number of corrector steps. 


(a) 


Momentum Predictor: Equation (106) is implicitly solved using the density 
and pressure from the previous time step. 


O 0,0 


V a p\ o * * o P y 

At + p° PUi = Ia nb U i,nb- VA iV + — ^T" 


(109) 


The resulting solution, u * , in general, will not satisfy continuity. 

(b) First Momentum Corrector: The momentum equation is written in an 

explicit corrector form as follows 


0 0 w° 

PU: V 


( V * ** ^ * vs A * > o 

— + — piq =Za nb u it 1 lb -VA i P + S U + —£ 


( 110 ) 


The momentum equation can be recast in delta form by subtracting Equation (109) 
from (110). 

rl 

. . ( * 

an) 


* ++ o * 


V a p 

P Ui-pUi=-\- 


+ 4 j VAiip'-p) 


However, we cannot solve for uf* because p*, and consequently p* are unknown. 
To find p*, we will enforce that uC* and p* satisfy continuity. That is 


A- (pX>- 


WOO 


P - P P_ 1 
At 'dtl v. 


( 112 ) 


* o o 


Taking a discrete divergence of Equation (111) and utilizing Equation (112), we get 

(113) 

— — V / 

Using the equation of state p* is evaluated as 


A {(| + J) *A } (p*-pT = A + 


Li£- + £- | 2 JL 

At At 1 v 
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p*=p* / (p“ T ‘) 


(114) 


substituting in Equation (113) 


A 


At 


+ 


v A 


A- (p°«;) 



(115) 


Solution to the above yields p*, which when substituted in Equation (111) yields the 
corrected velocity u.f*. 


(c) Enthalpy Predictor: The discretized enthalpy equation formulated with p*, p* 
and uC * is solved to yield h* or H*. From the known enthalpy, T* is obtained 
from thermodynamic property tables. 


(d) Second Momentum Corrector: The momentum equation is again 

represented in an explicit form 


p °u { V 
~At 


( \ y ^ p \ ** ** ** _ 

M + yj P =Ia nb U i,nb- VA iP +% + 

Recasting it in delta form by subtracting (110) from (116), 

P * - P U 7= (^ + ^j | £a nb[ u i*nb- u i,nb)- VA i (p" V)" a P 


f P- P 


V p 


U; 


(116) 


Enforcing continuity. 


p u i )=- 


** o Of yy 

P ~ P P 1 Vo 


At At 


1- 


(117) 

(118) 


Taking the discrete divergence of (117) and using (118) and the equation of state 


P 


/ ip". tO 


(119) 


we get the following equation for pressure correction 
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A V “p 

MU Vi 


1-1 


*4 it- 


f(p*,T*y 


At 


V - V 


1 


= 4 


t/ a p 


+ 


At p 


Sa n b{ U ~nb- “,‘j- a f ( P p f )“•“) + ^j{f (P*' T *)-f(p°' T )) 


( 120 ) 


Solving the above yields p** which can be used to evaluate p** from (119) and u*** 
from (117). 


The quantities u***, p**, p ** and T* are taken to be the flow variables at the new 
time level. Steps (a) through (d) are repeated to get the flow field for subsequent 
time steps. 


3.6 Crank-Nicholson Algorithm 


This algorithm is adopted for transient flow analysis and is formally second-order 
accurate in time. Using an appropriate weighting function, the algorithm can also 
be modified to give first-order accuracy in time. The second-order accuracy is 
achieved by evaluating all convective and diffusive fluxes as well as source terms at 
time level ( n+1/2 ) where n in the old time level. The algorithm consists of the 
following steps: 

a. Evaluate all flux and source terms at last time-step i.e. at time level n. 
These are stored till next time step. 

b. Intermediate velocity field, u*, v* are calculated using the momentum 
equations. The fluxes are calculated using the following expression 

(/c + /d) = a (/c + fot ^ ~ a X/c + /d)” 

where the superscript k denotes the iteration level, a is called the 
Crank-Nicholson parameter, and controls the time-implicitness and 
accuracy of the scheme. The method is second-order accurate for 
a =0.5. 
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c. Pressure corrections are evaluated as in the SIMPLEC algorithms, with 
the mass source term evaluated as 

m = o m k + (1 — a )m n 

The velocities and pressures are updated as in SIMPLEC. 

d. Steps b-c are repeated for a time step till a suitable convergence 
criterion is satisfied. All the flow variables are updated and 
calculations started for the next time step. 

By changing the value of the Crank-Nicholson parameter to 1.0, the method reduces 
to the Euler backward time-discretization which is 1st order accurate in time. 

3.7 Moving Grid Algorithm 

The moving grid algorithm is used to treat time-dependent flow problems with 
changing/ deforming grids as a result of boundary deformation, the formulation in 
SCISEAL is based on an extension of the Space Conservation Law (SCL) described by 
Demirdzic and Peric 1 * for 2-D problems. When the computational grid moves in 
time and space, the grid velocity, V g/ enters the flow equations and the relation 
between the grid velocity and the control volume is given as 

where V is the volume of control volume and S is its surface area. The continuity 
and the generalized transport equations for the problem are written in integral form 
as 

TtS> pdv+ i p ( v ' v *y ds=o (i22) 
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(123) 


d_ 

dt 



p<t{V-V g )d S = jq 



V 


where V is an arbitrary moving control volume, s is the surface of V and V is the 
absolute fluid velocity vector. In Equation (123), 0 could be the velocity vector or 
any other scalar quantity, q and S<p are the diffusive flux and source term for the 
corresponding flow quantity. The grid velocity Vg is an unknown, and the SCL is 
used to evaluate it as follows. 

Applying the space conservation law. Equation (121), to a moving control volume 
and using parametric time discretization, one obtains 


V- V> 

At 



V g ds + (l-0)j f 


V° g -ds 


o 


(124) 


where V stands for volume of the control volume at a new time step and 
superscript o stands for values at the old time step. Note that 0 = 1.0 is the standard 
Euler backward differencing and 6 - 0.5 is the Crank-Nicolson scheme and the value 
of 0 should be the same as the one used for other conservation equations for 
consistency. For numerical purposes, it is convenient to define a volumetric grid 
flux at a control volume (or cell) face, /, as 

(125) 

and therefore, the discretized SCL can be rewritten as 

(i 26 > 

If AVf stands for the volume swept by the cell face /, during one time step, it yields 
from the basic geometric requirement that 
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(127) 


V- V° 
At 


f 


f 


Hence, a natural way for calculating the volumetric grid flux is to let 


AV f 

-^ = ^ + ( 1 - 6 )^ ( 128 ) 

or 

AV f f i\ 

*/ = (ST'en (129) 

The apparent mass flux due to grid motion is obtained as 

Gj = P f & f (130) 

which when subtracted from the cell face mass flux, Gf, results in the effective mass 
flux that should be used in the evaluation of convective fluxes. The evaluation of 
Gf has been discussed in Section 3.4. 

A close look at A Vf reveals that AVf is simply the volume (but it could be positive or 
negative) and the conventional method of calculating the volume of a control 
volume is applicable also to A Vf. 

3.8 Domain Interface Treatment 

The multi-domain approach involves division of the overall flow domain into 
several suitable subdomains or zones such that all the grid generation in each of the 
subdomains is easier. Proper exchange of information across the interfaces where 
the subdomains join is crucial for the success of the solution procedure. The non- 
linear flow equations require a fully implicit and conservative treatment of the 
interface data exchange for an efficient and robust solution methodology. At present 
the SCISEAL code interface treatment allows a one-to-one match between cells as 
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well as a multi-to-one match, the only restriction being that an integer number of 
cells on one side of the interface must connect to a single cell on the other side. This 
is done through a special scheme in which the cell arrangement at an interface can 
be represented by a standard stencil shown in Figure 9. The bigger cell is identified 
as the parent (or base) cell while cells (one or many) on the other side of the 
interface are identified as children (or secondary cells). Such an arrangement allows 
the treatment of all domain interfaces in a standard form. For ease of illustration, 
the interface treatment is explained for 2D in the following sections. Extension to 
3D is fairly straightforward. 


S 2 


sA- 

Parent Cell “ 


— — • p + 


O p; ; 

•f 1 •Pi 

• 

O 

° A « 

f 2 •P2 

• 

h *P3 

• 

O 

U •pa 

• 

B 

• P* 



Child Cells 


Figure 9. Basic Interface Treatment Stencil 


3.8.1 Interface Interpolation 

In order to discretize the convection terms of the generalized convection-diffusion 
equation, the flow variable <p needs to be interpolated at cell faces. Consider the cell 
face at the domain interface of Figure 9. To evaluate the interface value at point fa 

(denoted by x), for example, we need to identify a fictitious point p 1 . If the <j> value 

at p 1 is obtained, the value at fa could be easily evaluated by any of the schemes 

described in Section 3.4. The p t value could be obtained in the following three 
ways: 
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(a) Upwind Method: 


(131) 



(b) 


Standard 2nd-Order: 


^ v\ 


Si 

& l + s 1 


<j>p + 


* I 


+ 


0p+ 


(132) 


(c) Flux-Limited Scheme: 


$p\ - Op + $p + ^ (^) 

f S p+ ) 

S p = S p _ Minimode 1,-p — 

l b P~! 

r 

Sp - = —faT 
c _<h>+-<i>p 
bp +~ Sp 


(133) 


Scheme (b) is the default for all variables except convective fluxes; the scheme for 
convective fluxes depends on the selection by user. 


3.8.2 Cross-Diffusion Term Treatment 

Cross diffusion terms are important when severely skewed BFC grids are used. 
They can be the dominant terms in flows with little or no convection as in heat 
conduction problems. In SCISEAL, the cross diffusion terms at the domain 
interfaces are treated implicitly to ensure accuracy. This treatment is illustrated for 
one term below 


^=7 IKS) (134) 

Integration of this term for the child cells pi, pz, ps and p 4 is essentially the same as 
that for a single-domain grid. The values of 0 at the child cell comers are evaluated 
from the surrounding nodal values, both physical and fictitious. For instance, 0^, 
the value of 0 at the northwest corner of cell p 4 , is estimated as 
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( 135 ) 


Integration of Equation (134) for the parent cell P requires special attention. 

[ [ D?J dt; dr! d(= (/r, 21 ^) - (/r* 21 ^) (136) 

J J Jp v ,J e v "w 

The second term on the right hand side of the preceding equation, i.e., the west face 
term is evaluated as explained in Section 3.3.5. The first term on the right hand side 
is evaluated at each of the 4 segments that form the east face of cell P. 



For example, for the segment fa, which is the interface between cell P and cell P 4 , 



(138) 


< p A and (ps are evaluated using relations such as Equation (135). 



In a number of seal and turbine secondary flow path problems the flow domain will 
contain internal solid objects, e.g. cavity wall, labyrinth seal teeth etc . Thermal 
energy transport can occur across the solid-fluid interfaces and through the solid 
objects. The energy equation solutions must consider convection-conduction in 
fluids, conduction in solids and the fluid-solid interface transfer. Such a procedure 
is called conjugate heat transfer (CHT) analysis. 


For thermal energy conservation, the following criteria have to be met at a solid- 
fluid interface (see Figure 10): 
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• Heat flux at the interface must be equal both on solid and fluid sides; 
and 

• Temperature has to be continuous across the interface. 

These two conditions can be mathematically represented as 


K s (VT) s -n = K f (V!) r n 

(139) 

ii 

5 ? 

h 

n 

b? 

(140) 


where the subscripts s and / stand for solid and fluid respectively, and the subscript i 
indicates the interface. Differences in the solid and fluid properties imply that 
quantities such as enthalpy will not be continuous across the interface. 



Figure 10. Interfacial Conditions for Conjugate Heat Transfer 
3.9.2 Equivalent Thermal Conductivity 

An equivalent thermal conductivity can be calculated using the interfacial balance 
equation. Equation (139) can be rewritten as 
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K 8 VT 8 -n = K f Vr r n = K e Vr.n 


( 141 ) 


where k e is the equivalent thermal conductivity. In tensor notation, yj . is 
expressed as 


VT-n = 


BT_ 

% 


3T e 2 oT e 
e 2 ■ e, + % e 3 


■ tt 


( 142 ) 


where e%, ei, and es are the covariant base vectors in 77 , and f directions, while e 1 , 
e 2 , and e 3 are the contravariant base vectors. 


By assuming 


ar s _T s -T,. 3T S T a -T fl 3T S T a - T a 
^ ' 3r? 4 t7 s ' 3f 


in the solid phase, and 


a Tf T r T f 3 T f T a - T a 3 T f Tq-T^ 

3£ ~ Afy ' dr, - Arj f flnd 3f " 4$ 

in the liquid phase, the temperature at the interface T* can be obtained and 
expressed as follows 


AT f +BT s + (D-C) + (F-E) 

X+B 


( 143 ) 


where 
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A = 


C = 


.ifi 


^ V b r n Jf 

H T n - T *) { - 2 


At} 


e 2 e 2 


if 


E = 


*7p3 " Tfc) I e~ • n 


Arf 


e 3 • e 


3j/ 


1C, 


B = 


D = 


Ui •« 


K s (T a -T 2 ) 


^ e 2 -n 


At] 


e 2 - 




F= 


M T s“ T *) 


( j „ 
1 e ■ n 


AC 


\e -e 3j 


Note that T^, T#, T& and T a are the interfacial temperatures at the corners. By 
substituting Equation (142) into Equation (141), the equivalent thermal conductivity 
K e is obtained after some manipulation. 


AB (T - T f + A(D +F) + B(C + E)) 
K * = (A + B)[G-(T s -r / + H + J)] 


where 


(144) 


G = 


A C s + A Cf V«J n 


,H = 


Tn-Ta 

( e 2 -it 

At] 

e 2 


, and l - 


Tu,-T 4 >- 3 




e • 7i 


e 3 • 


3.9.3 Turbulent Flow Considerations 

When turbulent flow with wall functions is being treated, the fluid conductivity kf 
needs to be modified to include the effects of turbulence. This is done by 
considering the expression for wall heat flux in turbulent flows 


--KVT • /i = 


p/ c *M T r T i) 

r 


(145) 


Here, T + is evaluated from the thermal law of the wall, and the friction velocity, u T , 
is 


u 


X 



(146) 
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An effective fluid conductivity K ^ is next defined as 


Pfpf U r( T f~ T i) 

T+ 


= K f VT‘n w 


(14 7 ) 


to yield the expression for Kj 


y 8 


where 8 is the normal distance between the wall and the first cell. 


( 148 ) 


Kf is used in place of in all the equations presented in Section 3.9. 
3.10 Turbulence Models 


A number of turbulence models have been incorporated in the SCISEAL code for 
the computations of the turbulent flows that exist in the seal flows. The theoretical 
framework behind turbulence modeling is outlined in this section, followed by a 
description of the various models that are available. The choice of the particular 
model that should be used in computations often will depend on factors such as the 
type of the flow, constraints on grid sizes, and flow details and any particular 
physical characteristics. 

3.10.1 Eddy Viscosity 

Favre averaging, the basis for all the turbulence models in SCISEAL, has already 
been described in Section 3.2.4. We have also seen how Favre averaging introduces 
additional terms known as Reynolds stresses in the Favre- Averaged Navier-Stokes 
(FANS) equations and how these stresses are modeled using the Boussinesq eddy 
viscosity concept (see Equations 27 and 28). Following the kinetic theory of gases, 
the eddy viscosity is generally modeled as the product of a velocity scale q and a 
length scale t 
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\i t = Cpq( 


(149) 


where C is a constant of proportionality. Various models differ in the way q and i 
are estimated and each of the following sections describes a turbulence model. In 
the description of models, the overbar for p and p, and tilde for u, v, etc. will be 
dropped for convenience. 

3.10.2 Baldwin-Lomax Model 

This belongs to the class of algebraic turbulence models because the velocity and 
length scales are obtained from algebraic relations. It is also commonly referred to as 
a mixing-length model because it employs Prandtl's mixing-length hypothesis in 
modeling length and velocity scales. 

Baldwin and Lomax 12 developed this model primarily for wall-bounded flows. Like 
the mixing-length model of Cebed and Smith 9 , it employs different expressions for 
fi t in the inner and outer parts of the boundary layer. 

_ / M-t inner V — V crossover 
\ M-f outer fo* — y crossover 


(150) 


In the inner layer, Prandtl's mixing-length model and the Van Driest's damping 
function are used to estimate the length scale 

Q = k y [1 - exp {-y+ / A + )\ (151) 

where A + is the Van Driest's damping constant. y + , the distance from the wall in 
wall units, is defined as 

y + = yu x /v, u z = jtjp (152) 

In the preceding expression u T is commonly known as the friction velodty with r w 
being the shear stress at the wall. 
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The velocity scale q is modeled as the product of p and the root mean square 
vorticity 


col = 


V9u _ M + ffe . + (to . Ml 

^dy dxj [dz dy j vox dz) 


1/2 


(153) 


Using the preceding expressions, the eddy viscosity in the inner layer is obtained as 


M't inner ~ P ^ 


(154) 


The outer layer eddy viscosity is determined from the following expression 


m outer ^ cp P ^ ’wake F kleb (]/) 


(155) 


where K is the Clauser constant, is an additional constant, and 

t r 


F wake ~ Ctlitl Yy^^x ^max > ^~ivk Vmax p 


Udif 


maxj 


(156) 


The quantities y max and F max are determined from the function 


%) = yM[l-«p(-y + M + )] 


(157) 


The quantity F max is the maximum value of F(y) that occurs within the boundary- 
layer and y max is the value of y at which the maximum occurs. 

F kleb (y) is the Klebanoff intermittency factor given by 


F]deb <y) 


1+55 


( Ckletif } 6 ] 1 

Umax J J 


(158) 


The quantity U di j is the difference between the maximum and minimum total 
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velocity in the boundary layer ( i.e at a fixed x station) 


u dif = vV + c 2 + a’ 2 ) mai - vV + tf 2 + u? 2 ) mj „ 


(159) 


The second term in tf^is zero for stationary walls. 

The values used for the constants appearing in the preceding expressions are 


A + =26 C =1.6 C kleb = 0.3 C wk = 0.25 k=0.4 K = 0.0168 


3.10.3 Standard k-e Model 

Several versions of k-e models are in use today, but the model employed in the seals 
code is based on Launder and Spalding 13 . This model employs two partial 
differential equations to estimate the velocity and length scales and hence 
commonly known as the two-equation model. These equations are the k- and e- 
equations which govern the transport of the turbulent kinetic energy (TKE) and its 
dissipation rate respectively. The modeled equations are 


^(pfc) + ^-(p^) = pP-pe + ^ 


V + Vt dk' 
&k dxj 


dt + dr • ( p U ) e ) ” Ce i P k ' C 


p e' 


3 x i 


^ + jte' 

a z dXj 


(160) 


with the production P defined as 


P = 



+ 


du i 9«; 2 t du 

1 T i -3d^ S 'ij-3F i -3 k lE 


m 

m 


(161) 


The square root of k is taken to be the velocity scale while the length scale is 
obtained from 
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cf it 3 ' 2 


The expression for eddy viscosity is 




~ T“ 


( 162 ) 


(163) 


The five constants used in the model are 

= 0.09 ; C Ej -1.44 ; C z =132; a k = 1.0; a E = 2.3 

The model uses wall functions to compute the turbulence quantities as well as the 
turbulent viscosity of the fluid at the computational points next to the walls. The 
wall functions described by Launder and Spalding 13 , are derived from experimental 
and analytical knowledge of the one-dimensional Couette flow which exists near 
the wall. A semi-empirical universal function of non-dimensional distance normal 
to the wall, y + , is 


y + 


p By- 

P 


(164) 


In the above definition, Sy is the distance normal to the wall and u T is the "friction 
velocity" given by 




i Vi 


(165) 


In the internal sublayer (y + > 11.63) the velocity variation may be described by a 
logarithmic relationship i.e. 


u^ln (Ey + ) 


(166) 


where E = 9.70 and k = 0.4034 are experimentally determined constants. 
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In both the viscous ( y + < 11.63) and internal (y + > 11.63) sublayers, the shear stress is 
calculated from the product of effective viscosity fi e ^ and normal velocity gradient 

dw/(9y, i.e. 


T 


w 



(1 67) 


where 


J p. for y + <11.63 
^ = U turb {ory + > 11.63 


(168) 


Near the wall, the transport equation for the turbulent kinetic energy, k, reduces to a 
balance between the local production and dissipation of k to give 



(169) 


The velocity gradient may be replaced. from Equation (168) and the dissipation rate 
from 


m = c,i p 


k 2 

x 


to give 


= c^ 2 pk = p u\ 


Hence, it follows from Equation (165) 


p Cf k 1 ' 2 a 
TW= y* { Ey+) 


(170) 


(171) 


(172) 
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Evaluation of the heat flux at the wall in the enthalpy equation is also done by 
assuming a profile for temperature near the wall. For laminar flows, the wall heat 
flux is obtained from finite differences because linear variation of temperature is a 
good approximation close to the wall. 



(173) 


For turbulent flows, Reynolds analogy is assumed, i.e. the velocity and temperature 
profiles are assumed to be similar when plotted in wall coordinates. The wall heat 
flux is obtained from the following expression 

C (T-T) + ^aJv 2 -V 2 ) 

J + _ pv w) 2 eff\ w) (174) 

q jp\ 


V and V w are the total velocities at the first grid-point and the wall respectively. T + 
at the first grid-point is obtained from a thermal law of the wall given below. 

T + = gu + fory + <y+ (175) 

= a { (w + P + ) for y + > y+ (176) 

In the preceding relations, P + is a parameter which is dependent on the fluid Prandtl 
numbers while y £ is the value of y + at which both the expressions (175) and (176) 
yield the same value of T. 



If o = 0.7 and a t = 0.9 we get P + - -2.13 and y* = 9.585. 


(177) 
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3.10.4 Multiple-Scale Model 

Most turbulence models, including the standard k-e model, assume one time scale 
for both the production and the dissipation rates of turbulence. However, experi- 
ments and the Direct Numerical Simulation (DNS) of turbulence have shown that 
most of the turbulence production occurs at large scales (energy carrying eddies) and 
is cascaded to smaller scales (dissipative eddies) where most of the dissipation takes 
place. Several authors ( e.g ., Hanjalic et a/., 14 ; Fabris et a/., 15 ) have suggested 
partitioning of the energy spectrum into production range k and dissipation range 

k t , the sum of the two quantities being k. Kim and Chen 16 have developed a 

multiple time-scale turbulence model based on variable partitioning of turbulent 
kinetic energy. The partitioning is determined as a part of the solution and depends 
on local turbulence intensity, production, energy transfer and dissipation rate. The 
partition is moved into higher wave numbers when production is high, and to low 
wave numbers when production vanishes. The model uses a transport equation for 
each k p/ k t , e p and e t , however these equations differ from each other only in the 

source terms. 


i ( pk ?) + ( p u > k v) = pp - p£ v + w; 

+ P“A) = P £ P -P E < + ^:(^S:' 


a pp 2 , pP£ p 

dt 


P e p 


d [ V + Pt dE V 


_a 

dt 


( p£ p) + dxj ( pU ^ ~ kp kp ' Ce P3 kp ° £ P dxj 

, A a , ^ ^ peF 2 P Pe p e f pe? , a (P + Pt 

( pe ^ + dx f lp u fr) ~ kt kt " C ^ kt ax ; [ ere t 'SxjJ 


(178) 


The eddy viscosity is defined as 


v t = C H k2 lep 


(179) 


where k, the total TKE, is calculated as 
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k = k +k t 


(180) 


The model constants are 

= 0.09 ) — Oft — 0.75 ; cr^ = c £ f = 1.15 

C^a =0.21; C £p2 =1.24; = 1 .84; C ell =0.29; C„ 2 =1.28; = 2.66 

3.10.5 Low Reynolds Number Model 

All three models described above, viz. the Standard k-e model, the extended k-e 
model, and the multiple time scale model, are of high-Reynolds number form and 
therefore require the used of wall functions. However, the commonly used wall 
functions may not be accurate in flows with large separation, suction, blowing, heat 
transfer, relaminarization, etc. This difficulty associated with wall functions can be 
reduced by the use of low-Reynolds number k-e models which permit integration of 
momentum and k-e equations all the way to the wall. The k-e equations are 
modified as shown below to include the effect of molecular viscosity in the near 
wall regions. 


^(pfc) + ^(pu y fc) = ^ 


7 


+ dk' 

dX; 
11 




o. {p£)+ o {pu . E)= 0 


3e 


+ C £l h 


+ p (P - e - D) 




(181) 


Several versions low-Reynolds number k-e models are available today. After a 
careful analysis of the comparative studies done by Patel et al} 7 , Brankovic and 
Stowers 18 , Avva et a/ 19 , the low-Re model of Chien 20 was selected for 
implementation in the seals code. The model parameters appearing in the 
preceding equations are: 

C^O.09; C Ei -135 ; C £i = 1.8 ; cr k = 1.0 ; <t e = 13 
f iL = l-exp (-0.0115 y + ) , f 2 = 1.0 ; f 2 = l - 0.22 exp [1RJ6) 2 ] 

D = 2vk/y 2 ; E = -2 v(e/y 2 ) exp (-0.5 y + ) 
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3.10.6 2-Laver Model 

The usual choice of turbulence models: the high- and low-Reynolds number models 
can be applied depending on the flow conditions and the grid sizes available. The 
model applicability depends on the distance of the first computational node from 
the wall. For the High-Re model, this distance, expressed in terms of the nondimen- 
sional distance y+, must be greater than 11.5 for the model to be accurate. On the 
other hand, the low-Re model required that the near-wall point be in the laminar 
sub-layer, at y+ of below 1. Computationally the high-Re model is robust and is the 
preferred choice. The low-Re model, under some circumstances becomes extremely 
stiff, and is harder to use. 

In seal flows, the narrow clearances typically can generate very low near-wall 
distances, even for relatively coarse grids. In such cases, when y+ falls below 11.5, 
the more robust standard k-e loses accuracy, and tends to overpredict the wall-shear. 
Use of the low-Re model in these cases will be more appropriate, however, the 
solutions are grid sensitive, and computationally stiff. In some cases, such as flows 
in labyrinth and grooved seals, there are regions where the flow is fast ( e.g . laby seal 
tip gap), while in some regions it is quite slow (laby cavity). In a flow of this type, 
even a coarse grid will generate very low y+ values in the cavity region, while the 
tip-gap region has relatively high near-wall y+ distances. In such a case, neither of 
the above models can be applied either due to loss of accuracy (high-Re model), or 
prohibitively large grid sizes (low-Re models). To treat such type of problems, a 2- 
layer model of wall treatment has been incorporated. The model provides the 
accuracy of the low-Re model for very low y+, and smoothly blends to a high-Re 
solution at higher y+ values. 

In the outer layer, the standard k-e equations govern the transport or turbulence and 
the eddy viscosity is computed as 

k 2 

I*r c ii7 (183) 


In the inner layer where the molecular viscosity is either dominant or comparable 
to eddy viscosity. Equation (183) is replaced by 
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m ~ c ^ \/k 


(184) 


The viscosity length scale l is determined using a Van Driest type correlation 

i» = Qy[ 1 - ex p{-¥)] 

Re is a local Reynolds number based on TKE 



(185) 


(186) 


For conformity with the log law, Cj is taken as 


r - * 

' “ c 3 / 4 


(187) 


where k is the von Karman constant. 

In the inner layer, the e-equation is replaced by an algebraic length scale equation 


e = 



(188) 


where the dissipation length scale 1 £ is modeled as 


/ = 

£ 1 + b/Re 


(189) 


Rodi 21 recommends the following values for constants a and b appearing in 
Equations (185) and (189): 


a = 50.5 


b - 5.3 
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To be consistent with the standard wall functions, k in Equation (187) is taken to be 
0.433. Note that a value of 0.40334 is used for k appearing the log law u + -1 /k (Ey + ). 


Equation (184) may be recast by using Equation (188) 


v t - */k ^ • e 

_ r i Jk ^ _ r l ± fc 2 _ r * fc 2 

'-ji. l \i ~TT 1 ~ v -'^ / ~ '“n M IT 

l £ l Z 


(190) 


Equation (190) resembles the typical formulation employed by most low-Re models. 
f is the ratio of viscosity length scale to dissipation length scale. Two-layer model is 
based on the premise that the viscosity length scale is smaller than the dissipation 
length scale inside the inner layer, but equal in the outer layer. 



1 - exp 



[ 1 

\i + 5.3/Re 


(191) 


Note that is strictly a function of Re only. The location where becomes unity is 
used in the present formulation as a matching criterion. 


Fie 

Numerator 
of Eq. (2.116) 

Denominator 
of Eq. (2.116) 

V 

1 

0.0196 

0.1587 

0.12348 

10 

0.1796 

0.6536 

0.275 

100 

0.862 

0.9497 

0.9076 

180 

0.97168 

0.9714 

1 .0003 

200 

0.981 

0.9742 

1.007 

1000 

0.9999997 

0.99473 

1 .0053 
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3.11 Solution Methods 


SCISEAL uses a pressure-based, iterative segregated solution procedure. The 
equations for the flow variables are solved in a sequential manner, and respectedly 
till a convergent solution is obtained. The overall solution method for SIMPLEC 
and PISO algorithm is presented in this section. 

3.11.1 SIMPLEC Algorithm 

The overall solution procedure for the SIMPLEC algorithm is shown in Figure 11. 
Note that all the parameters that dictate how many times a procedure is repeated 
can be specified by the user. These are the number of iterations (NITER), the 
number oc continuity iterations (CJTER) and, in the case of transient simulation, 
the number of time steps (NTSTEP). The procedure for a steady-state simulation is 
simply repeated at each time step of a transient simulation. The number of 
iterations to be performed is dictated by the overall residual reduction obtained. At 
each iteration the program will calculate a residual for each variable which is the 
sum of the absolute value of the residual for that variable at teach control cell. The 
residuals are not normalized, and hence the convergence criterion should be based 
on the reduction of residuals rather than on the absolute values. A drop of 4-5 
orders of magnitude is usually sufficient to consider the solution as converged. 
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Figure 11. Solution Flowchart for SIMPLEC Algorithm 
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3.11.2 PISQ Algorithm 

The overall solution procedure for the PISO algorithm is shown in Figure 12. 
Although PISO is an inherently transient procedure, it can be applied to steady state 
problems. In this case, the number of iterations (NITER) is comparable to the 
number of time steps performed in a transient simulation. The time increment 
through which the solution is allowed to advance at each time step is controlled by 
the under-relaxation factor applied and the time-asymptotic solution is sought. 



Figure 12. Solution Flowchart for PISO Algorithm 
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3.11.3 Under-Relaxation 

Under-relaxation of dependent and auxiliary variables is used to constrain the 
change in the variable from iteration to iteration in order to prevent divergence of 
the solution procedure. For all dependent variable except the pressure correction, 
this is done by modifying Equation (66) in the following way 

a p( I+/ )^ = - 2l ni>^i> +S U + (192) 

Where is the last iteration value. At convergence, and (p p are same and 
effect of I disapperars. During the iterative procedure, the factor I increases the 
diagonal term associated with <p p and in effect reduces the corrections made in <j>* p to 
get <j> . A larger value of / thus implies a higher relaxation. 

If we express I in the form of a transient term link coefficient 


/ = 


pV 

M f tt v 


(193) 


Then Atf can be viewed as a pseudo time step. In SCISEAL the value of I is specified 
directly and values between 0.2 and 0.8 are common. Higher relaxation may be 
needed for difficult problems. 

The auxiliary variables p, P, T and v can be under-relaxed by specifying a linear 
under-relaxation factor, a, which is applied in the following way 


er°= ae+{i- a)e 094 ) 

where 8 new is the updated value of the auxiliary variable, 8 is the value of that 
auxiliary variable that would be calculated with no under-relaxation and 8* is the 
current iteration value of the auxiliary variable. 

3.11.4 Linear Equation Solvers 

The discretized flow equation for each grid cell are interlinked with neighboring 
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cells, and when assembled generate a set of coupled linear equations. These can be 
written in the matrix form as 


m=s 

Due to the nature of the structured grid, the matrix A is sparse and has a diagonal 
banded structure. The individual non-zero elements of A are the link coefficients 
associated with each flow cell. The matrix A has to be inverted to get the solution 
vector 0. This can be done in SCISEAL using two types of linear equation solvers. 

3.11.4.1 Whole Field Solver Whole field solvers are based on the strongly implicit 
procedure (SIP). The method involves modification of the original linear equation 

set as 


[A + B]0 =S+[B]<j> (195) 

The matrix B is such that A + B is easily factored into upper and lower triangular 
matrices L and U. Further more, L and U matrices together have the same number 
of non-zero diagonals as the original matrix A. The modified equation then 

becomes. 


[L] [Id*- S + [B](f> 


This can also be expressed as 


[L]V = S + [B]0 


where 


V =[U]0 


(196) 


(197) 


(198) 


An iterative procedure can then defined and used 
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Evaluate V k + 1 from [L]V k + 1 =S + [B](f) k 
Evaluate <p k + 1 from [U\<p k + J = V k + 1 

The whole-field solvers in SCISEAL use this procedure. The form of L and U is 
such that points along lines in one of the coordinate direction are linked implicitly 
through tri-diagonal arrays. There are, therefore, three of the solvers named the 
WHOLE-X, WHOLE-Y and WHOLE-Z solvers. 

3.11.4.2 Conjugate Gradient Squared Solver Conjugate gradient type solvers have 
many advantages over classic iterative methods such as suitability for vectorization 
and the lack of user-specified parameters. The CGS algorithm has been incorporated 
into the SCISEAL program. The CGS algorithm applied to the system A<p = S is 
expressed as follows 

Initialization (n=0) r 0 = S - A<p 0 

<lo = V-i = 0; p.j = I 
T Pn 

Iteration {n SO) p n = r 0 r n , $ n = — 

f'n - 1 

u n~ r n + 4n 
Pn~ u n + Pn (*?n @n Pn- 1) 

v n = A Pn 

T Pn 

0 n ~ r o v n' a n ~ 

n 

tfn+l ~ u n~ a n v n 

r n + \ = r n - a n A { u n + % + l) 

3 = ^x -KXjfn + ln + l) 


The convergence rate of conjugate gradient algorithms depends on the spectral 
radius of the coefficient matrix and can be effectively accelerated by preconditioning 
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the system. This preconditioning is accomplished by transforming the system 
A<p = p'A(f> = p S . Incomplete Cholesky decomposition is used as a preconditioner 
in SCISEAL. 


3.12 



As a rotor rotates in a fluid seal, the fluid forces developed in the seal play an 
important part in the overall stability of the rotating machinery. The forces 
developed in a seal can be large enough to destabilize a rotor, or stabilize an 
otherwise unstable system. The seal force characteristics also play an important part 
in determining the critical frequencies of a rotating system. 


The seal flow forces generated when a rotor moves in the seal are related to the 
rotor center displacement, velocity and accelerations (Figure 13) through the 
following relation 



K 

yy 

-K 

*y 


K yz m , Cyy Cyz 
K L Z J -C C 

zz zy zz 



M M 

yy y 2 

-M M 

zy zz 



(199) 


where F y and F z are the fluid reaction forces as a result of rotor center motion, 
Kyy,K zz the direct stiffness coefficients, the cross-coupled stiffness coefficients, 

C yy ,C zz the direct damping coefficients, C yz ,C zy the cross-coupled damping 
coefficients, M yy ,M zz the direct inertia (mass) coefficients, and M yz ,M zy the cross- 
coupled inertia coefficients. These coefficients relate the reaction forces to the rotor 
center displacents, velocities and accelerations y, z, y z, y, z . The values and signs of 
these coefficients determine the effect of the fluid flow in the seal on the dynamic 
characteristics of the supported rotor system. 
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Figure 13. Rotor Displacement, Velocity and Acceleration and Fluid Reaction 
Forces in a Generic Seal Configuration 

At present there are two different methods available in the SCISEAL code to 
evaluate these coefficients. These are: 

a. whirling rotor method; and 

b. small-perturbation method. 

The whirling rotor method is based on computations that involve the full CFD 
solution of the seal flow . The perturbation method is based on computations of the 
perturbed fluid flow under a prescribed, small rotor motion. A description of these 
methods follow. 

3.12.1 Whirling Rotor Method 

In this method, the seal rotor is assumed to undergo a whirl such that the rotor 
center describes a circle about the stator center (Figure 14). Such a problem is time- 
dependent, but the circular rotor whirl orbit allows a coordinate-frame 
transformation that renders the problem quasi-steady. This is done by using a 
reference frame that has the origin at the stator center and is rotating with the rotor 
at its whirl speed. In this frame, rotor appears stationary, and a steady-state analysis 
is possible with the inclusion of appropriate body forces in the momentum 
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equations. A description of this transformation is given after the treatment for the 
perturbation method, as well as in Reference 45. 



Figure 14. Seal Configuration for the Circular Whirl Orbit Method for 
Rotordynamics 

Since the rotor whirls around the stator center, the system can generate only a skew- 
symmetric set of rotordynamic coefficients described as 



K k 
-k K 




M m 
-m M 



( 200 ) 


where K,k now are the direct and cross-coupled stiffness coefficients, C,c are the 
direct and cross-coupled damping coefficients, and M,m are the direct and cross- 
coupled inertia coefficients. These are assumed to be independent of the whirl fre- 
quency Q . Next, we use the definitions of the displacement of the rotor center 


y = r^os (at) 
z = r^in ( Q.t ) 


( 201 ) 


to calculate the rotor center velocities and accelerations, with r 0 as the radius of the 
whirl orbit. These are then used in Equation (200) to yield relations between 
reaction forces and various rotordynamic coefficients 
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( 202 ) 


F v . 

- / = K + cQ - MQ 2 

- = -k + CQ + mO 2 
7 0 

The calculation procedure then involves computations of the CFD solutions in the 
seal with whirling rotor at several different whirl frequencies (at least 3). The 
pressure fields on the rotor surface are integrated at each whirl frequency to provide 
F y and F z as functions of the whirl frequency. An appropriate curve-fit is used on 
the force vs. whirl frequency curves to yield the rotordynamic coefficients. 

As mentioned earlier, this method is applicable to the centered rotor position. 
When the rotor whirls about a non-centered position, even the use of a rotating 
frame of reference cannot make the problem quasi-steady. For this reason, the skew- 
symmetric dynamic coefficient set can be calculated using this method. 

3.12.2 Small Pert urbation Method 

This method is based on the solutions of perturbations in the seal flow field under a 
prescribed small motion of the rotor center about a given nominal position (similar 
to the method described by Nordmann and Dietzen 22 ). The perturbation variables 
can be solved using a quasi-steady solution method to yield the fluid reaction forces. 
The method is sufficiently general to include statically eccentric and/ or misaligned 
seals. 

In this method, the rotor of the seal is assumed to undergo a circular whirl around a 
given nominal position, with a very small orbit radius r 0 , given by 

r 0 = eC 0 (203) 

where C 0 is a reference length, e.g. the nominal seal clearance, and e is a small 
number. In the nominal position, the rotor can be centered, eccentric and/or 
misaligned. A time-dependent flow field exists in the seal as a result of the rotor 
motion, and the flow variables are assumed to have the form, e.g., 
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u = u 0 + eu : 
v = v + evj 


( 204 ) 


where u Q , v 0 correspond to the steady-state seal flow solutions with the rotor in the 
nominal position (0th order), and u v i are the time-dependent perturbation 
solutions (1st order). Using these definitions in the Navier-Stokes equations, 
ignoring terms with e 2 and higher terms, and separating out the 0th and 1st order 
terms, the governing equations for the two sets of variables are generated. Since the 
rotor center motion is a combination of sine and cosine functions of time, the 1st 
order flow variables are also assumed to have the form, e.g., 

u i =u ic cos (Of) + u ls sin (Of) 

v i = v ic cos (Q*) + v is sin (Q*J (205) 


where u lc , u ls etc. are functions of space only. Similar definitions for all flow 

variables are assumed and substituted in the 1st order flow equations. By separating 
out the terms containing sin(£2t) and cos(£2t) in each of the equations, two equations 
for each of the 1st order flow quantity can be obtained ( i.e . one equation each for 
quantities like u lc , u ls and so on). 

The next step involves regrouping and redefining new 1st order variables such as: 


uj = u u + i u ls 

v t = u lc + i v ls (206) 

where i =v c T . The equations for these complex quantities are generated in a similar 
fashion: 


(Equation for Uj) - ( Equation for u lc ) + i (Equation for u ls ) (207) 

The resulting equations have a form very similar to the 0th order equations. Thus, 
the momentum equation for has the form: 
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( 208 ) 


-iQ (p„ u 3 + p 3 w 0 ) / + A (/Po V* • e l ) + iJPo “ Ao • e l ) 

dq os 

+ (/ P:“o • et ) +7 Vt • e k ) = ^ (w* |t) + S ul 

where S ul contains the complex source terms that arise due to the complex 

convective and diffusive fluxes as well as the grid transformations. The continuity 
and energy equations are transformed in a similar fashion, and both contain 
additional complex source terms. Several points to note here are: a) The 

momentum equations are linear, since the convective terms contain the Oth order 
fluxes which are constant; 2) the equations are similar in form to the base Oth order 
equations, so that a similar solution procedure can be used and 3) the 1st order 
variables are complex quantities, so that complex algebra is required for their solu- 
tions. 

The perturbations are assumed to be small enough so that the turbulence quantities: 
k, e and the turbulent viscosities are assumed unchanged in the 1st order equations. 
This reduces the 1st order set to only the momentum, continuity and energy equa- 
tions (for compressible flows) only. 

The overall solution procedure is: 

a. calculation of the steady-state, Oth order solutions of the Navier-Stokes 
equations for the seal flow with the rotor in the nominal position; 

b. solution of the 1st order equations using the Oth order convective 
fluxes and turbulent quantities; and 

c. use the resulting pressure field to calculate the rotordynamic 
coefficients. 

Once the 1st order pressure field is known, it is integrated over the surface of the 
rotor and the resulting time dependent forces are decomposed into Y and Z 
components to yield F yr , F yi , F zr , and F zi as given below: 

Fyr = [ Vic ■ njrdedx) (209a) 
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Fyi = U | Pu ■ ”y( rdQdx) 

(209b) 

F Z x = _[ Pu ■ (rdQdx) 

(209c) 

F zr = ^| | Pic ' n^rdQdx) 

(209d) 


where n y and n z are the components of the local surface normals in the y and z 
Cartesian directions. Using the definition of the displacement of the rotor center 


y -r 0 cos (ftf) 
z = r 0 sin (Qf) 


( 210 ) 


the rotor center velocity and accelerations are calculated. These are then substituted 
together with the definitions in Equations (209) in the basic force-displacement 
relation. Equation (199). Collecting terms with cos(Qt) and sin(Qt), four equations 
linking the integrated reaction forces and the whirl speed Q. are generated, as given 

below: 

~F yr = Kyy - Cy Z Cl - MyyO. 

-F yi = K yz -c yy n-M yz n 2 
-F zr = -K^ + C 22 Q - M zy il 2 

— F ZI - = Xjj. + C Z yl2 ~ iVl zz ii (211) 

As in the whirling rotor method, the solutions of the 1st order quantities are 
generated at several different whirl frequencies Q, and appropriate curves are fitted 
to the F yr , Fy i etc. vs. frequency curves to yield all the rotordynamic coefficient. 


As stated earlier, this method generated the full, non-symmetric set of coefficients, 
and can be extended to generate the coefficients due to angular displacements. The 
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overall solution procedure requires the solutions of several quasi-steady problems, 
and hence is computationally fast. 

3.12.3 Grid Transformation for Concentric Whirling Rotor 

When a rotor whirls inside a seal, the moving surface of the rotor generates a 
deforming grid which produces a time-dependent flow. This flow can be solved 
using the moving grid method described earlier. This problem is time-dependent 
and hence computationally expensive. In a special case of rotor whirl, however, it is 
possible to change the frame of reference and render the problem quasi-steady. This 
special case refers to the rotor whirl which has a circular orbit and stator center as 
the orbit center. This transformation is the basis of the whirling rotor method and 
following is a brief description of the transformation and the resulting flow 
problem. 

The whirling rotor is shown in Figure 15. The spin speed is co x rad/s and the whirl 
speed is Q x rad/s. The figure shows the rotor when its center is in orbit at one time 
instant. The orbit radius is e and r T is the radius vector of a typical point on the rotor 
surface with respect to the rotor center. 

To render this problem quasi-steady, one must switch the reference frame such that 
in the transformed frame, the flow domain appears non-deforming even when it 
deforms in the absolute reference frame. When a rotor is whirling in a circular 
orbit, this transformed frame corresponds to a frame that is rotating at rotor whirl 
speed, has the stator center as the origin, and stator axis as the axis of rotation. 

To understand this transformation, refer to Figure 16. In this figure the positions of 
the rotor center at two different time instances are shown. Also shown are two local 
reference frames such that the local Y' axis always joins the stator center and the 
minimum clearance point. In the local Y' - Z' frame the flow domain now appears 
unchanged, and hence Y' - Z' is the reference frame we need. Since the minimum 
clearance point moves at whirl speed Q x , so does the Y' - Z' plane. 
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Figure 15. Schematic of a Spinning Rotor Whirling in a Circular Orbit about the 
State- a rnter 



Figure 16. Positions of the Whirling Rotor at Two Time Instants. Also Shown are 
the 'Local' Reference Frames / - z' at these Two Instances. 
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The changes in the flow equations as a result of this rotating frame transformation 
involve addition of the centrifugal and Coriolis accelerations to the momentum 
equations. The wall boundary conditions also need to be transformed. This needs 
care and the procedure is outlined below. Referring to Figure 17, the absolute wall 
velocities are: 

stator V - 0 (212) 

Sfd 

and 

rotor V =(3 x/ 1 + 5 xe (213) 

ta x r x 

where f , § etc. refer to the position vectors of surface points. For the rotor surface 

the total velocity is a combination of the rotor spin and whirl speeds. 
Transformation to the rotating frame involves subtraction of a velocity of the type 

V t =Q xR (214) 

0 X 

from all wall velocities, where R refers to the radius vector with respect to the stator 
center as shown in Figure 17. With this transformation the surface velocities 
relative to the transformed frame become: 


Stator: 

V -V -Q xR = -Q xR 

rjr s,a x s x s 

(215) 

Rotor: 

V ~ V -Q x R = c5 x f + Q x2-Q x R 

r/ rA x r x r x x r 

(216) 

Noting that 

R = f + e 

t r 

(217) 

we get the rotor velocity as V f - (co x - Q J x ? r 

(218) 


Note that in the transformed frame, the stator wall appears to move backwards. The 
rotor wall, which describes a whirl in the absolute frame, has been reduced to a 
spinning wall such that the modified spin speed is the difference between the rotor 
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spin and whirl speeds. This corroborates our assertion that in the rotating frame the 
grid deformations are removed and the problems can be treated as quasi-steady. 



Figure 17. Definitions of the Rotor Spin and Whirl Velocities as well as Radius 
Vectors for Wall Velocity Transformation. Note that all Definitions 
are in the absolute Reference Frame. 

Interpretation of the results obtained using this frame must be done with care. 
With the transformation the eccentric rotor looks similar to a bearing, and there is a 
tendency to explain the flow in terms of the rotor and stator wall absolute velocities. 
This procedure will lead to unphysical interpretations. The proper way to analyze 
this flow is either to treat it as a time-dependent flow (absolute reference frame) or 
use the transformed plane, with proper transformations to the wall velocities, 
where the stator wall has a motion. 


NAS A/CR— 2004-2 1 3 1 99/VOL6 


81 



4.0 SAMPLE AND VALIDATION PROBLEMS 


A large number of flow problems have been simulated with SCISEAL and the 
results compared with benchmark analytical /experimental /numerical date to assess 
and validate the accuracy of the numerical and physical models incorporated in the 
code. A partial list of the problems is given below. Following the list, descriptions 
of several relevant problems and selected results are given. The first several 
problems in the list are checkout type 2D problems and hence details are omitted. 
All seal-related and other relevant problems are described. 

4.1 Problem Titles 

1. Fully-developed flow in a pipe and channel. 

2. Developing laminar flow in a narrow annulus between two cylinders. Slug 
flow at inlet, fully-developed flow at outlet. 

3. Laminar flow between rotating cylinders. Below critical Taylor number, 
tangential flow only. 

4. Flow between two cylinders, rotating inner cylinder. Taylor vortex flow. 
Laminar and turbulent 23 . 

5. 2-D driven cavity flow, Reynolds number up to 10,000. Comparisons with 
numerical results by Ghia et.al , 24 . 

6. Couette flow under different pressure gradients. With and without heat 
transfer. 

7. Planar wedge flow in a slider bearing. 

8. Laminar flow over a back step. Reattachment length comparison with 
experiments by Armaly and Durst 25 . 

9. Shock reflection over a flat plate. 

10. Turbulent flow in a plane channel. Fully-developed solution at exit 
compared with experiments by Laufer 26 . 

11. Turbulent flow induced by rotating disk in a cavity. Comparison with 
experiments by Daily and Nece 27 . 

12. Centripetal flow in a stator-rotor configuration. Comparison with 
experiments by Dibelius et.al. 28 . 

13. Flow between stator and whirling rotor of a seal. 2-D results for 0, 0.5, and 
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synchronous whirl frequencies. 

14. Flow over a bank of tubes. 

15. Turbulent flow in an annular seal. Comparison with experiments by 
Morrison et.al . 29 . 

16. Turbulent flow in a 7-cavity labyrinth seal. Comparison with experiments by 
Morrison et.al. 30 . 

17. Turbulent compressible flow and heat transfer in turbine disk cavities 31 . 

18. Laminar flow in a square duct with a 90° bend. Comparison with 
experimental data by Taylor et.al. 31 . 

19. 3-D driven cavity flow with lid clearance and axial pressure gradient. Control 
of flow through vortex imposition 33 . 

20. Flow in cavities on a rotor for an electrical motor. Interaction of Taylor 
vortices with driven cavity flow. 

21. Flow in infinite and finite length bearings (without cavitation) 

Comparison of calculated attitude angles with theory 34 ' 35 . 

22. Flow and rotordynamic coefficient calculation for straight, incompressible 
seals. Comparison with results from other numerical and analytical 

solutions 36 . 

23. Flow and rotordynamic coefficients in tapered compressible flow seals. 
Comparison with bulk-flow theory results 37 . 

24. Rotordynamic Coefficients in a long annular incompressible flow seal. 
Comparison with experimental data 38 . 

25. Calculation of entrance loss coefficients in the entrance region of a generic 
seal. Effect of flow and geometry on the loss coefficient values 39 . 

26. Flow coefficient and pressures in a 5 cavity, straight knife, look-through 
labyrinth seal. Comparison with experimental data 40 . 

27. Flow coefficients and pressures in a 2 cavity, tapered knife, look-through 
labyrinth seal. Comparison with experimental data 41 ' 42 . 

28. Flow coefficients and pressures in a 2 cavity, straight-knife, stepped labyrinth 
seal. Comparison with experimental data 41 ' 42 . 

29. Computations of the rotordynamic coefficients for an eccentric annular seal 

with incompressible flow 43 . 

30. Flow solutions in a whirling annular seal 44 and comparison with 
experiments. 45 " 47 
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31. Numerical solutions of flow and cooling effectiveness in rim seals and disc 
cavities 48 and comparison with experiments. 49 

32. Interaction of mainpath-secondary flow in multiply connected turbine 50 disc 
cavities and comparison with experimental data. 51 

33. Flow and conjugate heat transfer simulations in turbine disc cavities. 52 
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4.2 Problem Description and Sample Results 


Following are the descriptions and SCISEAL results and comparisons with other 
date /results where applicable. The number associated with each of the problems 
refers to it's position in the list given above. 

11. Turbulent Flow Induced bv Rotating Disk in a Cavity 
Problem specification 

• Calculation of the flow induced by a rotating disk in an enclosed cavity. 
Benchmark data 

• Experimental measurements from Daily and Nece 27 . 

Grid 

• 40 cells in the axial direction, 60 cells in the radial direction with 
clustering near the walls. 

Boundary conditions 

• Specified angular velocity for the rotor walls. 

• Wall conditions for all other boundaries. 

Numerics and physical models 

• Central differencing with 0.05 damping. 

• Standard k-e model with wall functions. 

Results 

• Flow geometry as shown in Figure 18a. 

• Normalized radial and tangential velocities at a given radius are 
shown in Figures 18b and 18c. Also shown in the figures are the 
experimental data from Daily and Nece 27 . 
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□ Data of Daily & Neca 
— Colocated Grid 


x/b 



(b) Radial Velocity (c) Tangential Velocity 

Figure 18. Turbulent Flow Due to a Rotor in an Enclosed Cavity. Experimental 
Data from Daily and Nece 27 
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15. Annular Seal Flow 
Problem specification 

• Calculation of turbulent flow in an annular seal. 

Experimental data 

• Experimental data by Morrison, et al. 29 . 

Grid 

• 25 cells in the radial direction, 58 cells in the axial direction; cells in 
radial direction clustered near the walls. 

Boundary conditions 

• Experimental profiles of the velocities and turbulence quantities at 
inlet boundary. 

• Specified pressure at the outflow boundary. 

• Wall condition with specified angular speed at rotor wall. 

• Stationary wall conditions at stator wall. 

Numerics and physical models 

• Central differencing with 0.01 damping. 

• Standard two equation k-e model for turbulence. 

Results 

• Geometry of the rotor is shown in Figure 19a, and the experimental 
setup is shown in Figure 19b. 

• Computed and experimental contours of the axial, azimuthal and 
radial velocities are shown in Figures 20, 21 and 22, respectively. 

• Figure 23 shows the computed turbulent kinetic energy profiles. 
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Figure 19. 




94 0 mm 


17 J mm 


(a) Geometry of Annular Rotor 



(b) Annular Seal Inlet Geometry 

Flow Details for the Annular Seal. Reynolds number based on the gap 
= 27000, Taylor number = 6600, shaft speed = 3600 rpm 
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Figure 20. Annular Seal Flow. Contours of the scaled axial velocity u x /u. (a) near 
inlet (0 < x/c < 4), (b) near outlet (25 < x/c < 29.4) 
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Figure 21. Annular Seal Flow. Contours of the scaled azimuthal velocity, 
u 0 / w shaft (a) near inlet (0 < x/c <4), (b) near outlet (25 < x/c < 29.4) 
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Experimental Numerical 

Figure 22. Annular Seal Flow. Contours of the scaled radial velocity u r /u. Near 
inlet details (0 < x/c < 4) 



(b) 


Figure 23. Annular Seal Flow. Contours of the scaled turbulent kinetic energy, 
j (' u r + u q Numerical results, (a) near inlet (0 < x/c < 4), (b) 

near outlet (25 < x/c < 29.4). 
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16. Seven Cavity Labyrinth Seal 

Problem specification 

• Calculation of turbulent flow in a seven-cavity labyrinth seal. 
Experimental data 

• Experimental data by Morrison, et al 30 . 

Grid 

• 30 cells in the axial and radial directions per cavity, 

• 10 cells in the radial clearance between the rotor tooth and the stator. 

• Stretching used to cluster the grid near the rotor and stator walls. 

Boundary conditions 

• Experimental profiles for velocities and turbulence quantities at inlet 
boundary. 

• Specified pressure at outflow boundary. 

• Wall condition with specified angular velocity at rotor walls. 

• Wall conditions at stator wall. 

Numerics and physical models 

• Central differencing with 0.01 damping. 

• Standard two equation k-e model for turbulence. 

Results 

• Details of the rotor are shown in Figure 24a, and the experimental 
setup is shown in Figure 24b. 

• Computed and numerical velocity vector plots are shown in Figure 25. 

• Computed and experimental contours of the axial, radial and 
tangential velocities are shown in Figures 26, 27 and 28, respectively. 

• Figure 29 shows computed contours of the turbulent kinetic energy. 
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3.048 mm 



(a) Details of the Labyrinth Rotor 



(b) Seal Inlet Geometry 

Figure 24. Seven Cavity Labyrinth Seal Flow. Reynolds number based on the 
clearance = 28000, taylor number = 7000, shaft speed = 3600 rpm. 
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Experimental 


Numerical 


Third Cavity 

Figure 26. Seven Cavity Labyrinth Seal. Contours of scaled axial velocity, u x /U. 
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First Cavity 



Experimental Numerical 

Third Cavity 

Figure 27. Seven Cavity Labyrinth Seal. Contours of scaled radial velocity, u r /u. 
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First Cavity 



Experimental Numerical 

Third Cavity 

Figure 28. Seven Cavity Labyrinth Seal. Contours of scaled azimuthal velocities, 
V w 5 haft- 
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First Cavity 


Third Cavity 


Figure 29. 


Seven Cavity Labyrinth Seal Flow. Contours of normalized turbulent 
7 ( '2 '2 ' 2 \ / 

kinetic energy 2 l M r + u q +u x)/ if- Numerical results. 
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17. Turbulent Compressible Flow in Turbine Disk Cavities 
Problem Specification 

Computation of the turbulent compressible flow is representative turbine disk 
cavities, under the effect of secondary cooling flows. The effect of varying cooling 
flow rates and geometry changes are assessed. 

Geometry, Grid, Physical Models and Results 

Please refer to Reference 31 for detailed description of the problem and results. 
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19. 3-D Driven Cavity Flow with Lid Clearance and Axial Pressure Gradient 

Problem Specification 

To compute the interaction of the tip vortex and the core flow in a turbine blade 
passage. The problem is simulated using a driven cavity flow with lid clearance and 
axial pressure gradient. Effects of a vortex imposed upstream of the flow are 
assessed for vortex strength, sense and location. 

Grid, Geometry, Physical Models and Results 

Please refer to Reference 31 for a detailed description of the problem and results. 
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20. Flow in Cavities on a Rotor in an Electrical Motor 
Problem Specification 

Computation of the flow in the channels present in an electrical motor, to assess the 
flow structure, generation of Taylor vortices and their interaction with the flow in 
the channels. 

Grid, Geometry, Physical Models and Results 

Please refer to Reference 33 for a detailed description of the problem and results. 
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Computation of flow in A) an infinite length bearing, and B) finite length bearing 
without cavitation. 

Benchmark Data 

Analytical results for the bearing attitude angle as given in Fuller 34 , (infinite length) 
and Cameron 35 (finite length). 

Grid 

A. Infinite length case 

• Rotor radius = 36 mm, nominal clearance = 0.036 mm. 

• 40 cells along the circumference, 10 cells in the radial direction, 3 cells 
in the axial direction. 

B. Finite length case 

• Rotor radius = 36 mm, nominal clearance = 0.036 mm, axial length = 36 
mm (L/D = 0.5). 

• 10 cells in the axial direction, 40 along the circumference, 10 cells in the 
radial direction. 

Conditions 

Symmetry conditions on boundaries in axial direction. 

Rotating wall on rotor, (speed = 5,000 rpm), stationery wall on stator. 

Specified ambient pressure on both boundaries in axial direction. 
Rotating wall on rotor, (speed = 5,000 rpm), stationery wall on stator. 



NAS A/CR— 2004-2 1 3 1 99/VOL6 


101 



Physical models 

• Laminar, incompressible flow. 

• Central differencing with 0.1 damping. 

Results 

• Results for bearing eccentricities from 0.1 to 0.95 obtained for both cases 
A and B. 

• Half Sommerfeld conditions used to calculate rotor pressure forces and 
the resulting attitude angles. Attitude angles compared with analytical 
results. 

• Figure 30 shows attitude angles for Case A, together with analytical 
results from Fuller. 34 Figure 31 shows results for Case B, with 
analytical results from Cameron 35 , very good /excellent agreement seen 
in both cases. 
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Figure 30. Attitude Angles for Case A, Long Bearing 
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22. Rotordynamic Coefficients for Incompressible Annular Seal 

Problem Specification 

Computation of the rotordynamic coefficients for a short annular grid with 

incompressible flow. 

Experimental Data 

Experimental and numerical data from Dietzen and Nordmann. 36 

Grid 

• Seal radius = 23.5 mm, seal clearance = 0.2 mm, seal length = 25 mm. 

• 10 cells in axial direction, 30 in circumference direction and 5 or 12 in 
seal gap depending on the turbulence model used. 

Boundary Conditions 

• Specified upstream /downstream pressure differential with entrance 
loss coefficient = 0.5. 

• Stationary wall on stator, rotating wall on rotor, rotor speed varied 
from 1000-5000 rpm. 

Physical Models 

• Central-differencing with 0.1 damping. 

• Data obtained with standard k-e model with wall functions and the 2- 
layer model for turbulence. 

• Whirling rotor, numerical shaker and perturbation model tried for 
rotordynamic coefficients. 


NAS A/CR— 2004-2 1 3 1 99/VOL6 


104 



Results 

• Rotordynamic coefficients calculated at several rotor speeds. 

• Numerical results for whirling rotor method shown here. These 
results match very well with results from the numerical shaker 
method and perturbation method. 

• Figures 32 and 33 show the direct and cross-coupled stiffness 
coefficients. Figures 34 and 35 show the direct and cross-coupled 
damping coefficients and Figure 36 shows the direct inertia coefficient. 
Also shown are results from experiments, finite difference methods 
and bulk-flow theory. Very good agreement is observed in the present 
results and other published results. 
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Nordmann, f.d. 

- — — k-e model with wall f. 

2-layer model, y+ 0{1 ) 

— 2-layer model, y+ 0(20) 



Shaft speed, rpm. 


Figure 34. Direct Damping Coefficient, Annular Seal 
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Figure 35. Cross Coupled Damping Coefficient, Annular Seal 
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23. Rotordvnamic Coefficients for Tapered Seal with Compressible Flow 
Problem Description 

Computation of compressible flow and the rotordynamic coefficients for a tapered 
annular seal. 

Data 

Bulk flow theory results by Nelson. 37 
Grid 

• Seal radius = 32.5 mm, inlet and exit clearance = 0.172 and 0.086 mm 
L/D ratio varied from 0.1 to 0.4. 

• 12 cells in axial, 6 cells in radial and 30 cells in circumferential 
direction. 

Boundary Conditions 

• Upstream pressure = 1.52 MPa, exit pressure = 0.65 MPa, y = 1.4. 

• Upstream Temperature = 650k and entrant loss factor ^ as given in 
Reference 37. 

• Rotating wall on rotor stationery wall on stator, shaft speed = 30,400 
rpm. 

Physical Models 

• Central differencing with 0.1 damping. 

• Standard k-e model with wall functions. 

• "Numerical Shaker"' method for rotordynamic coefficient calculations. 
Results 

• Results obtained at L/D = 0.1, 0.2 and 0.4. 

• Results shown in Table 2.1 together with bulk flow theory results from 
Nelson 37 and good agreement between the two sets is seen. 
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Table 2. Tapered Gas Seal Rotordynamic Coefficients 


L/D 

KN/m 

k N/m 

C N-s/m 

c N-s/m 

Exit Mach 
Number 

0.1 

• ' ; . 948000 


\ 113 

' 0*016 


1150000 

15429 

9.94 

0.043 

1.00 

0.2 

1700000 

» 75700 

43*9 

' \ * : 0.073 

jg 0*96 

2125500 

60886 

38.62 

0.09 

0.97 

0.4 


" 267000 


; / .0*27 

0.83 

3553200 

233820 

145.9 

0.57 

0.83 


Nelson, 1985 37 


Present Results 
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24. Rotordynamic Coefficients of a Long Annular Seal 

Problem Specification 

Computation, of rotordynamic coefficients for a long annular seal with 

incompressible flow. 

Experimental Data 

• Experimental results reported by Kanemori and Iwatsubo. 38 

Grid 

• Seal radius = 39.656 mm, clearance = 0.394 mm, length = 240 mm. 

• 18 cells in axial direction, 30 in circumferential direction, 16 or 18 cells 
in seal gap depending on the axial pressure differentials. 

Boundary Conditions 

• Specified upstream /downstream pressure differential ranging from 20 
to 900 KPa with entrance loss factor \ - 0.5. 

• Rotating wall at rotor surface, stationary wall on stator; shaft speed = 
600, 1080, 1980 and 3,000 rpm. 

Physical Models 

• Central differencing scheme with 0.1 damping. 

• 2 turbulence models used (1) low-Re k-e (ii) 2 layer. 

• Problem tried with (i) whirling rotor method and (ii) small 
perturbation method for rotordynamics. 

Results 

• Similar results with both turbulence models and rotordynamics 
calculation methods. 

• Results from low Re k-e model and whirling rotor method reported 
here. 

• Figures 37 and 38 show the direct and cross-coupled stiffness 
coefficients. Figures 39 and 40 the direct and cross-coupled damping 
and Figure 41 the direct inertia (mass) coefficient. 

• Experimental data also plotted for comparison, very good agreement 
between computational and experimental 38 results. 
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Direct Stiffness Coefficient 



"Figure 37. Direct Stiffness Coefficients, Long Annular Sea? 


Cross-Coupled Stiffness Coefficient 



Figure 38. Cross Coupled Stiffness Coefficients, Long Annular Seal 
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Direct Damping Coefficient 



Figure 39. Direct Damping Coefficient, Long Annular Seal 
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Figure 40. Cross Coupled Damping Coefficient, Long Annular Seal 
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Figure 41. Direct Mass Coefficient, Long Annular Seal 
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25. Calculation of Entrance Loss Coefficients 


Entrance loss factor is an empirical factor used in seal calculations to account for the 
flow losses as the fluid enters the narrow seal clearance from a relatively large gap in 
the entrance region. A generic seal entrance region is considered and several 
geometry and flow parameters are varied to assess their effect on the loss factor. 

Problem Specification 

Compute the entire loss factors for a generic seal-entrance region and assess effects of 
1) rotor radius- to-clearance ratio, 2) Entrance region radial width-to-clearance ratio 
and 3) Flow rate on the entrance loss factor. 

Grid 

• Seal radius 25 mm, seal length = 25 mm, inlet region length 1.5 times 
the seal radius. 

• 30 cells in entrance region, 20 cells in the seal in axial direction. 5 cells 
in the seal clearance radially, and 30 or 50 cells radially in entrance 
region, depending on entrance-gap-to-clearance ratio. 

Boundary Conditions 

• Fixed downstream pressure, fully-developed turbulent axial velocity 
profile at inlet of entrance region for specified mean axial velocity. 

• Rotating wall on the rotor surface, stationary wall on stator. 

• Inlet k & e fully-developed profiles. 

Physical Models 

• Incompressible flow. 

• Central-differencing with 0.1 damping. 

• Standard k-e model with wall functions. 

Results 

• Typical grid shown in Figure 42. 

• Rotor speed at 3,000 rpm, rotational Re = 2,500 - 7,500. 

• Entrance loss coefficients at different clearance-to-radius ratios, 

Reynolds meters and Entrance-gap-to-clearance ratios shown in 


NAS A/CR— 2004-2 1 3 1 99/VOL6 


115 



Table 3. 

• Numerical results show: 

Radius-to-clearance ratio has the highest effect ratio? => ^ T 
Gap-to-clearance ratio T => § T with lower sensitivity 
for a fixed geometry, axial Re T => § i , with high sensitivity 
£ range from 0.406 to 0.68. 


Entrance Seal 



Figure 42. Flow Domain and One of the Grids Used for the Entrance Loss 
Coefficient Calculations 
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Table 3. Entrance Loss Coefficients, Radius /Clearance = 50 


Entrance Gap/Clearance = 50 

Entrance Gap/Clearance = 100 

m/s 

Re ax 

4 

u ax m/s 

Re ax 

\ 

10.814 

10377 

0.471 

10.82 

10384 

0.490 

16.232 

15484 

0.431 

10.24 

15584 

0.488 

21.619 

20746 

0.414 

21 .66 

20785 

0.482 

26.942 

25854 

0.406 

27.06 

25970 

0.48 


PA-92-24 (2 


Entrance Gap/Clearance = 50 

Entrance Gap/Clearance = 100 

Uax tti/s 

Re ax 


Uax rtVs 

Re ax 


10.80 

5181 

0.562 

10.797 

5167 

0.567 

16.56 

7945 

0.54 

16.176 

7761 

0.558 

21.595 

10361 

0.526 

21.55 

10339 

0.55 

26.67 

12796 

0.51 

26.934 

12664 

0.54 

32.27 

15484 

0.493 

32.24 

15469 

0.537 

43.062 

20667 

0.478 

42.533 

20408 

0.524 
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Entrance Gap/Clearance = 100 

m/s 

^ax 

% 

Uax m/s 

Reax 


10.82 

3461 

0.66 

10.75 

3438 

0.68 

16.19 

5178 

0.65 

16.09 

5146 

0.66 

21.49 

6874 

0.647 

21.47 

6874 

0.65 

26.74 

8553 

0.637 

26.81 

8553 

0.648 

32.25 

10315 

0.628 

32.176 

10292 

0.64 

48.33 

15461 

0.606 

47.87 

15315 

0.63 

64.487 

20630 

0.595 

64.165 

20630 

0.624 


PA-92-24 14 
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26. Five Cavity Straight-Through Labyrinth Seal with Straight Knives 


Experimental data on mass flow coefficients at different seal clearances and seal 
pressure ratios are reported in Reference 40. Also reported are detailed pressure 
profiles along the seal length at one pressure ratio and different seal clearances. 
Computations were performed at three different pressure ratios and three clearance 
values for the flow coefficients comparison. Comparison of pressure profiles was 
done at four different clearance values. The computed results were compared with 
experimental data and excellent agreement was obtained. Flow geometry details and 
results of the calculations are discussed below. 

Problem Specification 

Computation of steady-state turbulent air flow in a 5 cavity, look-through, planar 
labyrinth seal at different pressure ratios and tip gaps. 

Experimental Data 

Experimental results by Willig et a/. 40 
Grid 

30 x 30 cells in each cavity, 8 to 12 cells in the tip gap, 30 cells in axial direction in the 
upstream and downstream, portions. Appropriate clustering of cells near walls. 

Boundary Conditions 

• Downstream pressure kept fixed, upstream pressure varied to match 
specified pressure ratio across the seal. 

• Stationary walls on stator and rotor. 

• Inlet k and e values specified using estimated inlet velocities. 

Physical Models 

• Compressible flow. 

• Central differencing scheme with 0.1 damping. 

• k-e model of turbulence with wall functions. 
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Results 

• Flow geometry and parameters are shown in Figure 43. 

• Flow solutions obtained at several pressure ratios across the seal, and at 
different tip gap values. 

• Figure 44 shows the pressure profiles in the seal for a pressure ratio 
p u /p d of 1.38 for four tip gaps. Computational and experimental values 

plotted for comparison. 

• Figure 45 shows the non-dimensional mass flow coefficient (<j>) (see 
Reference 40 for definition) as a function of the seal pressure ratio at 
different tip gaps. 



Figure 43. Details of the Flow Geometry for Five Cavity Look-Through Labyrinth 
Seal, b - 2.5, h = 10.5, t = 12, s = 0.5 - 2.52 (all dimensions in mm) 
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Figure 44. Non-Dimensionalized Pressure Drop at the Centers of the Cavities for 
a Pressure Ratio of 1.38 



Figure 45. Mass Flow Coefficient (<j>) as a Function of Pressure Ratio Across the 
Seal 
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27. Two Cavity, Straight-Through Labyrinth Seal with Tapered Knives 


Experimental 41 and computational 42 results for this case are available. 
Experimental data includes mass flow coefficients at different pressure ratios across 
the seal, and detailed pressure measurements in the seal for a fixed pressure ratio. 
Computations were performed at three different pressure ratios for comparison of 
flow coefficients as well as pressure profiles. Very good agreement was obtained 
between the computed and experimental pressure profiles, while the flow 
coefficient results show a good comparison, with some under prediction in 
computed results. Computational details and the results are given below. 

Problem Specification 

Computation of steady-state turbulent air flow in a two cavity, look-through, planar 
labyrinth seal with tapered knives. 

Experimental Data 

Experimental data by Tipton, et fl/. 41 

Grid 

• 26 cells in the axial and 50 cells in the radial direction per cavity and 210 
cells in the gap. 

• 51 and 75 cells in the axial direction in the upstream and downstream 
directions, respectively. 

• Grid is clustered near the rotor and stator walls. 

Boundary Conditions: 

• Exit pressure fixed, upstream pressure varied to match specified 
pressure ratio across the seal. 

• Inlet k and e values specified using estimated inlet velocities. 

• Wall conditions are specified along the stator and rotor walls. 

Numerical and Physical Models: 

• Central differencing with 0.1 damping. 

• k-e model of turbulence with standard wall functions. 
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Results: 


• Details of the flow domain are shown in Figure 46. 

• Figures 47 and 48 show the pressure drops along the direction of flow 
(as a function of inlet pressure) measured near the stator and rotor 
surfaces, respectively. The inlet/ exit pressure ratio is 2.0. 

• The leakage mass flow rates through the seals for different inlet /exit 
pressure ratios are plotted in the form of a flow parameter, <j> (see 
Reference 41 for definition) in Figure 49. 



Figure 46. Flow Domain for the 3-Knife Labyrinth Seal. All dimensions are in 
inches. Upstream and downstream region lengths are 3 and 5 inches, 
respectively. 
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Figure 47. Pressure Drop Along the Stator Surface Plotted as a Function of Inlet 
Pressure Along the Seal Length. A: mid-point of 1st knife; B: mid- 
point of 1st cavity; C: mid-point of 2nd knife; D: mid-point of 2nd 
cavity; E: mid-point of 3rd knife; F: half of the cavity width 
downstream of 3rd knife. Inlet /Exit pressure ration = 2.0 
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Figure 48. Pressure Drop Along the Rotor Surface Plotted as a Function of Inlet 
Pressure Along the Seal Length. A: mid-point of 1st knife; B: mid- 
point of 1st cavity; C: mid-point of 2nd knife; D: mid-point of 2nd 
cavity; E: mid-point of 3rd knife; F: half of the cavity width 
downstream of 3rd knife. Inlet /Exit pressure ration = 2.0 
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Figure 49. Leakage Mass Flow Rates 
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28. Two Cavity, Stepped Labyrinth Seal with Rectangular Knives 


As in test case 26, experimental and computational data is available on flow 
coefficients and pressure profile in the seal for the stepped seal. Computations were 
done at three different pressure ratios for flow coefficients as well as pressure values 
in the seal. The comparison of flow coefficients is good, with somewhat higher flow 
coefficients predicted by the CFD code. This could partly be due to the upstream and 
downstream boundary condition definitions. The comparison of the pressures in 
the seal is very good. Details of the computations are shown below. 

Problem Specification 

Computational of steady-state turbulent air flow in a two cavity, stepped labyrinth 
seal with rectangular knives. 

Experimental Data 
Experimental data by Tipton, et a/. 41 

Grid 

• 26 cells in the axial direction in each cavity, 24 radial cells before the 
step and 53 radial cells after the step in the first cavity; 43 radial cells 
before the step and 62 radial cells after the step in the second cavity. 

• 10 cells in the radial clearance and 5 cells in the axial direction between 
each knife tip and the stator wall. 

• A 40 x 44 grid in the upstream region and a 60 x 52 grid in the 
downstream region. 

• Grid is clustered near the rotor and stator walls. 

Boundary Conditions 

• Exit pressure fixed, upstream pressure varied to match specified ratio 
across the seal. 

• Inlet values of k and e are calculated using estimated inlet velocities. 

• Wall conditions are specified along the stator and rotor walls. 

Numerical and Physical Models 

• Central differencing with 0.1 damping. 
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k-e model for turbulence with standard wall functions. 


Results 




Figure 50. 


Details of the flow domain are shown in Figure 50 
Figures 51 and 52 show the static pressures along the direction of flow 
(as a function on inlet pressure) measured near the stator and rotor 
surfaces, respectively. The inlet/ exit pressure ratio is 2.0. 

The variation in flow parameter, <j> (see Reference 41 for definition), 
with the inlet/exit pressure ratio is plotted in Figure 53. 

O.l 

i 

n "T" 

[— 0 . 69 — | | 


0.125 


1.5 


2.1 


Details of the Flow Geometry of the Stepped Labyrinth Seal. All 
dimensions are in inches. Inlet and exit region lengths are 3 and 5 
inches, respectively. 
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Figure 51. Pressure Drop Along the Stator Surface Plotted as a Function of Inlet 
Pressure Along the Seal Length. A: location upstream of 1st knife 
equidistant to location C from location B; B: mid-point of 1st knife; C: 
mid-point of region in 1st cavity before step; D: mid-point of region in 
1st cavity after step; F: mid-pionot of region in 2nd cavity before step; 
G: mid-point region in 2nd cavity after step; H: mid-point of 3rd knife; 
I: location downstream of 3rd knife equidistant to location G from 
location H. Inlet /exit pressure ration = 2.0 
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Figure 52. Pressure Drop Along the Rotor Surface Plotted as a Function of Inlet 
Pressure Along the Seal Length. B: mid-point of 1st knife; C: mid-point 
region in 1st cavity before step; D: mid-point of region in 1st cavity after 
step; E: mid-point of 2nd knife; F: mid-point of region in 2nd cavity 
before step; G: mid-point of region in 2nd cavity after step; H: mid- 
point of 3rd knife. Inlet/ exit pressure ratio = 2.0 
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Figure 53. Flow Parameter as a Function of Inlet/Exit Pressure Ratio 
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29. Rotordvnamic Coefficients for an Eccentric Annular Seal 


Problem Specification 

Computations of rotordynamics coefficients for an annular seal and incompressible 

flow at different static eccentricities. 

Experimental Data 

Experimental and numerical results as reported by Simon and Frene. 43 

Grid 

• Seal radius = 80 mm, seal clearance = 0.36 mm, seal length = 40 mm. 

• 10 cells in axial direction, 5 in radial gap and 30 in circumferential 
direction. 

Boundary Conditions 

• Specified upstream and downstream pressure values with an extreme 
loss factor £ = 0.5. 

• Stationery wall on stator, rotating wall on rotor with shaft speed = 4,000 
rpm. 

Physical Models 

• Central differencing with 0.1 damping. 

• Standard k-e model with wall functions. 

• Perturbation method for rotordynamic coefficients. 

Results 

• Rotordynamics coefficients computer at static eccentricities ranging 
from 0 to 0.7. 

• Direct stiffness K ^ and K zz as well as K yz and K zy shown in Figures 54 

through 57. Also shown are experimental results and other numerical 
results for comparison. Very good agreement between present results 
and experimental results is seen. 
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• Damping coefficients Cyy, C^, and C zy shown in Figures 58 through 
61. Comparison with other numerical results also shown, very good 
agreement is seen between the data sets. 

• Inertia (mass) coefficient and shown in Figures 62 and 63. 
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Figure 54. Direct Stiffness, K , Annular Eccentric Seal 



Figure 55. Direct Stiffness, IC Z2 , Annular Eccentric Seal 
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Figure 56. Cross Coupled Stiffness, K^, Annular Eccentric Seal 



Figure 57. Cross Coupled Stiffness, K^, Annular Eccentric Seal 
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Figure 58. Direct Damping, Cyy, Annular Eccentric Seal 



Figure 59. Direct Damping, C^, Annular Eccentric Seal 
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Figure 61. Cross Coupled Damping, C zy , Annular Eccentric Seal 
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Figure 62. Direct Inertia, Myy, Annular Eccentric Seal 
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Figure 63. Direct Inertia, Annular Eccentric Seal 
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30. Flow in a Whirling Annular Seal 


Problem Specification 

Computation of the flow in a whirling annular seal to investigate the behavior and 
comparison with experimental velocity and pressure measurements. 

Grid, Geometry, Physical Models and Results 

Please refer to Reference 45 for detailed description of the problem, boundary 
conditions and results. 
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31. Flow and Cooling Effectiveness of Rim Seals 


Problem Specification 

Evaluation of the flow in a generic disc cavity and rim seal configuration. Use 
passive scalar transport of a tracer gas to evaluate the cooling effectiveness of four 
different types of rim seals at different flow conditions and compare with 
experiments. 

Grid, Geometry, Physical Models and Results 

Please refer to Reference 48 for a detailed description of the problem and results 
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32. Interaction of Secondary Flow and Mainpath Flow in Multiple Pi&c Cavities 


Problem Specification 

Analysis of the flow dynamics in a multi-cavity turbine disc configuration. Account 
for the coupling between cavities and main path flow. Use tracer gas transport to 
identify flow features and compare with experimental date. 

Grid, Geometry, Physical Models and Results 

Please refer to Reference 50 for a detailed description of the problem and results 
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33. Flow and Conjugate Heat Transfer Analysis in Turbine Disc Cavities 


Problem Description 

Simulation of flow and heat transfer in the turbine disc cavities of an actual turbine 
engine (Allison T-56). Evaluate mass flow rates and gas temperatures at various rim 
seals, effect of interstage labyrinth seal and conjugate heat transfer in labyrinth seal 
support discs. 

Grid. Geometry, Physical Models and Results 

Please refer to Reference 52 for a detailed description of the problem and results 
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